Note
Go to the end to download the full example code.
Benchmark Shi, Xie, Xuan & Nocedal's general method¶
The goal of this example is to benchmark the ShiXieXuanNocedalGeneral
class on a collection of test problems.
These problems are created by the build_benchmark()
static method, which returns a list of problems.
import numpy as np
import pylab as pl
import tabulate
import numericalderivative as nd
Compute the first derivative¶
class ShiXieXuanNocedalGeneralMethod:
def __init__(
self,
relative_precision,
initial_step,
differentiation_order=1,
formula_accuracy=2,
direction="central",
):
"""
Create a ShiXieXuanNocedalGeneral method to compute the approximate first derivative
Parameters
----------
relative_precision : float, > 0, optional
The relative precision of evaluation of f.
initial_step : float, > 0
The initial step in the algorithm.
"""
self.relative_precision = relative_precision
self.initial_step = initial_step
self.differentiation_order = differentiation_order
self.formula_accuracy = formula_accuracy
self.direction = direction
def compute_derivative(self, function, x):
"""
Compute the first derivative using ShiXieXuanNocedal
Parameters
----------
function : function
The function
x : float
The test point
Returns
-------
f_derivative_approx : float
The approximate value of the first derivative of the function at point x
number_of_function_evaluations : int
The number of function evaluations.
"""
formula = nd.GeneralFiniteDifference(
function,
x,
self.differentiation_order,
self.formula_accuracy,
self.direction,
)
absolute_precision = abs(function(x)) * self.relative_precision
algorithm = nd.ShiXieXuanNocedalGeneral(
formula,
absolute_precision,
)
step, _ = algorithm.find_step(self.initial_step)
f_derivative_approx = algorithm.compute_derivative(step)
number_of_function_evaluations = algorithm.get_number_of_function_evaluations()
return f_derivative_approx, number_of_function_evaluations
The next example computes the approximate derivative on the
ExponentialProblem.
initial_step = 1.0e0
problem = nd.ExponentialProblem()
print(problem)
function = problem.get_function()
x = problem.get_x()
differentiation_order = 1 # First derivative
formula_accuracy = 2 # Order 2
formula = nd.GeneralFiniteDifference(
function,
x,
differentiation_order,
formula_accuracy,
direction="central", # Central formula
)
algorithm = nd.ShiXieXuanNocedalGeneral(
formula,
verbose=True,
)
third_derivative = problem.get_third_derivative()
third_derivative_value = third_derivative(x)
optimal_step, absolute_error = nd.FirstDerivativeCentral.compute_step(
third_derivative_value
)
print("Exact h* = %.3e" % (optimal_step))
function = problem.get_function()
first_derivative = problem.get_first_derivative()
x = 1.0
relative_precision = 1.0e-15
absolute_precision = abs(function(x)) * relative_precision
method = ShiXieXuanNocedalGeneralMethod(absolute_precision, initial_step)
(
f_prime_approx,
number_of_function_evaluations,
) = method.compute_derivative(function, x)
first_derivative_value = first_derivative(x)
absolute_error = abs(f_prime_approx - first_derivative_value)
print(
"x = %.3f, error = %.3e, Func. eval. = %d"
% (x, absolute_error, number_of_function_evaluations)
)
DerivativeBenchmarkProblem
name = exp
x = 1.0
f(x) = 2.718281828459045
f'(x) = 2.718281828459045
f''(x) = 2.718281828459045
f^(3)(x) = 2.718281828459045
f^(4)(x) = 2.718281828459045
f^(5)(x) = 2.718281828459045
Exact h* = 4.797e-06
x = 1.000, error = 4.193e-10, Func. eval. = 63
Perform the benchmark¶
print("+ Benchmark on several points")
number_of_test_points = 21 # This number of test points is odd
problem = nd.SinProblem()
print(problem)
interval = problem.get_interval()
function = problem.get_function()
first_derivative = problem.get_first_derivative()
initial_step = 1.0e2
test_points = np.linspace(interval[0], interval[1], number_of_test_points)
relative_precision = 1.0e-15
absolute_precision = abs(function(x)) * relative_precision
method = ShiXieXuanNocedalGeneralMethod(absolute_precision, initial_step)
average_relative_error, average_feval, data = nd.benchmark_method(
function,
first_derivative,
test_points,
method.compute_derivative,
verbose=True,
)
print("Average relative error =", average_relative_error)
print("Average number of function evaluations =", average_feval)
tabulate.tabulate(data, headers=["x", "Rel. err.", "F. Eval."], tablefmt="html")
+ Benchmark on several points
DerivativeBenchmarkProblem
name = sin
x = 1.0
f(x) = 0.8414709848078965
f'(x) = 0.5403023058681398
f''(x) = -0.8414709848078965
f^(3)(x) = -0.5403023058681398
f^(4)(x) = 0.8414709848078965
f^(5)(x) = 0.5403023058681398
x = -3.142, abs. error = nan, rel. error = nan, Func. eval. = -1
x = -2.827, abs. error = 2.221e-11, rel. error = 2.335e-11, Func. eval. = 87
x = -2.513, abs. error = 3.156e-11, rel. error = 3.901e-11, Func. eval. = 93
x = -2.199, abs. error = 5.525e-11, rel. error = 9.399e-11, Func. eval. = 75
x = -1.885, abs. error = 2.995e-11, rel. error = 9.691e-11, Func. eval. = 75
x = -1.571, abs. error = 8.482e-17, rel. error = 1.385e+00, Func. eval. = 21
x = -1.257, abs. error = 2.762e-11, rel. error = 8.937e-11, Func. eval. = 75
x = -0.942, abs. error = 5.757e-11, rel. error = 9.795e-11, Func. eval. = 75
x = -0.628, abs. error = 3.814e-11, rel. error = 4.715e-11, Func. eval. = 93
x = -0.314, abs. error = 2.453e-11, rel. error = 2.580e-11, Func. eval. = 87
x = 0.000, abs. error = nan, rel. error = nan, Func. eval. = -1
x = 0.314, abs. error = 2.453e-11, rel. error = 2.580e-11, Func. eval. = 87
x = 0.628, abs. error = 3.814e-11, rel. error = 4.715e-11, Func. eval. = 93
x = 0.942, abs. error = 5.757e-11, rel. error = 9.795e-11, Func. eval. = 75
x = 1.257, abs. error = 2.762e-11, rel. error = 8.937e-11, Func. eval. = 75
x = 1.571, abs. error = 8.482e-17, rel. error = 1.385e+00, Func. eval. = 21
x = 1.885, abs. error = 2.995e-11, rel. error = 9.691e-11, Func. eval. = 75
x = 2.199, abs. error = 5.525e-11, rel. error = 9.399e-11, Func. eval. = 75
x = 2.513, abs. error = 3.156e-11, rel. error = 3.901e-11, Func. eval. = 93
x = 2.827, abs. error = 2.221e-11, rel. error = 2.335e-11, Func. eval. = 87
x = 3.142, abs. error = nan, rel. error = nan, Func. eval. = -1
Average relative error = nan
Average number of function evaluations = 64.71428571428571
Notice that the method does not perform correctly for the point \(x = 0\) for the polynomial problem. This test point appears only if the number of test points is odd, because the test interval is symmetric with respect to \(x = 0\). For this problem, the method does not perform correctly because the value of the function is zero at \(x = 0\). The method can perform correctly in this case, if it is provided a consistent value of the absolute error of the function value. Here, we compute the absolute error depending on the relative error and the absolute value of the value of the function. If the value of the function is zero, then the computed absolute error is zero, which produces a failure of the method.
Map from the problem name to the initial step.
initial_step_map = {
"polynomial": 1.0,
"inverse": 1.0e0,
"exp": 1.0e-1,
"log": 1.0e-3, # x > 0
"sqrt": 1.0e-3, # x > 0
"atan": 1.0e0,
"sin": 1.0e0,
"scaled exp": 1.0e5,
"GMSW": 1.0e0,
"SXXN1": 1.0e0,
"SXXN2": 1.0e0,
"SXXN3": 1.0e0,
"SXXN4": 1.0e0,
"Oliver1": 1.0e0,
"Oliver2": 1.0e0,
"Oliver3": 1.0e-3,
}
The next script evaluates a collection of benchmark problems
using the ShiXieXuanNocedalGeneral class.
number_of_test_points = 100 # This value can significantly change the results
data = []
function_list = nd.build_benchmark()
number_of_functions = len(function_list)
average_absolute_error_list = []
average_feval_list = []
relative_precision = 1.0e-15
delta_x = 1.0e-9
for i in range(number_of_functions):
problem = function_list[i]
name = problem.get_name()
initial_step = initial_step_map[name]
function = problem.get_function()
first_derivative = problem.get_first_derivative()
interval = problem.get_interval()
lower_x_bound, upper_x_bound = problem.get_interval()
print(f"Function #{i}, {name}")
if name == "sin":
# Change the lower and upper bound so that the points +/-pi
# are excluded (see below for details).
lower_x_bound += delta_x
upper_x_bound -= delta_x
test_points = np.linspace(lower_x_bound, upper_x_bound, number_of_test_points)
method = ShiXieXuanNocedalGeneralMethod(relative_precision, initial_step)
average_relative_error, average_feval, _ = nd.benchmark_method(
function, first_derivative, test_points, method.compute_derivative
)
average_absolute_error_list.append(average_relative_error)
average_feval_list.append(average_feval)
data.append(
(
name,
average_relative_error,
average_feval,
)
)
data.append(
[
"Average",
np.nanmean(average_absolute_error_list),
np.nanmean(average_feval_list),
]
)
tabulate.tabulate(
data,
headers=["Name", "Average rel. error", "Average func. eval"],
tablefmt="html",
)
Function #0, polynomial
Function #1, inverse
Function #2, exp
Function #3, log
Function #4, sqrt
Function #5, atan
Function #6, sin
Function #7, scaled exp
Function #8, GMSW
Function #9, SXXN1
Function #10, SXXN2
Function #11, SXXN3
Function #12, SXXN4
Function #13, Oliver1
Function #14, Oliver2
Function #15, Oliver3
Notice that the method cannot perform correctly for the sin function at the point \(x = \pm \pi\). Indeed, this function is such that \(f(x) = 0\). In this case, the absolute error is zero if it is computed from the relative error. Therefore, we make so that the points \(\pm \pi\) are excluded from the benchmark. The same problem appears at the point \(x = 0\). This point is not included in the test set if the number of points is even (e.g. with number_of_test_points = 100), but it might appear if the number of test points is odd.
Notice that the method cannot perform correctly for the polynomial problem because the third derivative is zero. This produces a zero test ratio which makes the method to fail.
Total running time of the script: (0 minutes 0.245 seconds)