set function
f
(
x
)
f(x)
f(x), at
[
a
,
b
]
[a,b]
[a,b] is second-order continuous differentiable and has a unique minimum point
x
0
x_0
x0?. because
f
(
x
)
f(x)
f(x) is
[
a
,
b
]
[a,b]
A unimodal function on [a,b], so
f
‘
‘
(
x
)
>
0
f”(x)>0
f”(x)>0,
x
∈
(
a
,
b
)
x\in(a,b)
x∈(a,b). to given
x
k
∈
[
a
,
b
]
x_k\in[a,b]
xk?∈[a,b], the function
f
(
x
)
f(x)
f(x) at
x
k
x_k
The nearest approximation of xk? is
q
k
(
x
)
=
f
(
x
k
)
+
f
‘
(
x
k
)
(
x
?
x
k
)
+
1
2
f
‘
‘
(
x
k
)
(
x
?
x
k
)
2
.
q_k(x)=f(x_k) + f'(x_k)(x-x_k) + \frac{1}{2}f”(x_k)(x-x_k)^2.
qk?(x)=f(xk?) + f′(xk?)(x?xk?) + 21?f′′(xk?)(x?xk?)2.
because
q
k
‘
(
x
k
)
=
f
‘
(
x
k
)
q’_k(x_k)=f'(x_k)
qk′?(xk?)=f′(xk?),
q
k
‘
‘
(
x
k
)
=
f
‘
‘
(
x
k
)
>
0
q”_k(x_k)=f”(x_k)>0
qk′′?(xk?)=f′′(xk?)>0, so
x
k
?
f
‘
(
x
k
)
f
‘
‘
(
x
k
)
x_k-\frac{f'(x_k)}{f”(x_k)}
xkf”(xk?)f'(xk?)? is
q
k
(
x
)
q_k(x)
The stagnation point of qk?(x) is also the only minimum point.
We use the following strategy to calculate
f
(
x
)
f(x)
The optimal solution of f(x)
x
0
x_0
Search sequence for x0?: take
x
1
∈
[
a
,
b
]
x_1\in[a,b]
x1?∈[a,b] is sufficiently close to
x
0
x_0
x0?, from
k
=
1
k=1
k=1 starts to iterate, if
f
‘
(
x
k
)
=
0
f'(x_k)=0
f′(xk?)=0, by
f
(
x
)
f(x)
f(x) at
[
a
,
b
]
[a,b]
Differentiability and unimodal knowledge on [a,b], find the optimal solution
x
0
=
x
k
x_0=x_k
x0?=xk?, as shown in (a) above. Otherwise, we use
q
k
(
x
)
q_k(x)
The unique minimum point of qk?(x)
x
k
?
f
‘
(
x
k
)
f
‘
‘
(
x
k
)
x_k-\frac{f'(x_k)}{f”(x_k)}
xkf”(xk?)f'(xk?)? as the first
k
+
1
k + 1
k + 1 iteration points
x
k
+
1
x_{k+1}
xk + 1?, ie
x
k
+
1
=
x
k
?
f
‘
(
x
k
)
f
‘
‘
(
x
k
)
.
x_{k + 1}=x_k-\frac{f'(x_k)}{f”(x_k)}.
xk + 1?=xkf”(xk?)f'(xk?)?.
no matter
x
k
x_k
xk? is in
f
(
x
)
f(x)
The rising range of f(x) is as shown in (b) above or in
f
(
x
)
f(x)
The decline interval of f(x) is shown in (c) above,
x
k
+
1
x_{k+1}
xk + 1? are expected to compare
x
k
x_k
xk? closer
f
(
x
)
f(x)
The minimum point of f(x)
x
0
x_0
x0?. Repeatedly, until the
f
‘
(
x
k
)
=
0
f'(x_k)=0
f'(xk?)=0.
The search algorithm based on this strategy is called Newton’s method, and the following code implements Newton’s algorithm.
from scipy.optimize import OptimizeResult
def newton(fun, x1, gtol, **options):
xk=x1
f1,f2=der_scalar(fun,xk)
k=1
while abs(f1)>=gtol:
k + =1
xk-=f1/f2
f1,f2=der_scalar(fun,xk)
bestx=xk
besty=fun(bestx)
return OptimizeResult(fun=besty, x=bestx, nit=k)
Lines 2-12 of the program define the function newton that implements Newton’s algorithm. The parameter fun represents the function
f
(
x
)
f(x)
f(x), x1 represents the initial point
x
1
x_1
x1?, gtol means error tolerance
ε
\varepsilon
ε, options are used to make minimize_scalar pass actual parameters such as x1 and gtol to newton. Line 3 initializes the iteration point xk to x1. Line 4 calls the derivative calculation function der_scalar for details, see the blog post “Numerical Calculation of the Derivative of a Unary Function”) calculation
f
(
x
)
f(x)
The first and second derivatives of f(x) are placed in f1 and f2. Line 5 sets the iteration count to 1. The while loop in lines 6 to 9 performs iterative calculation: the number of iterations k in line 7 is incremented by 1, and the iteration point is calculated in line 8
x
k
?
f
‘
(
x
k
)
f
‘
‘
(
x
k
)
x_k-\frac{f'(x_k)}{f”(x_k)}
xkf”(xk?)f'(xk?)? update xk. Line 9 calls der_scalar again to calculate the first and second derivatives at the updated iteration point. Repeatedly, until
∣
f
‘
(
x
k
)
∣
<
ε
|f'(x_k)|<\varepsilon
∣f′(xk?)∣<ε. Lines 10 and 11 respectively calculate the optimal solution approximation bestx and the optimal solution function approximation besty. Value line 12 returns the calculation result besty (optimal solution function value
f
(
x
0
)
f(x_0)
Approximate value of f(x0?)), bestx (optimal solution
x
0
x_0
x0?) and k (the number of iterations) to construct an OptimizeResult type object as the return value.
Example 1 Use the newton function defined by the above program to calculate the function
f
(
x
)
=
x
2
+
4
cos
?
x
f(x)=x^2 + 4\cos{x}
f(x)=x2 + 4cosx in
x
1
=
1.5
x_1=1.5
The minimum value point near x1?=1.5.
Solution: The following code completes the calculation.
import numpy as np #import numpy
from scipy.optimize import minimize_scalar #import minimize_scalar
f=lambda x:x**2 + 4*np.cos(x) #Set objective function
minimize_scalar(f, method=newton, options={<!-- -->'x1':1.5,'eps':1.48e-8}) #calculate the optimal solution
It is not difficult to understand this program by using the comment information in the code. run the program, output
fun: 2.316808419788213 nit: 7 x: 1.8954942647118507
Readers can compare the results of this operation with the results calculated by using the dichotomy method for the same function in Example 1 in the blog post “Unary Function Search Algorithm – Dichotomy”. At the same level of error tolerance, the efficiency of Newton’s method (using only 7 iterations) is obviously higher than that of bisection method (using 27 iterations).