1 Definitions related to numerical derivatives
In this section, we overview the key definitions that relate to the
functionality of numDeriv.
A derivative is the instantaneous rate of change of a function (assuming its differentiability): \[ f'(x) = \frac{\mathrm{d} f}{\mathrm{d} x} := \lim_{h\to0} \frac{f(x+h) - f(x)}{h} \]
A numerical derivative is an approximation of the expression above for a small fixed h. For a scalar function, this approximation can be calculated with two function evaluations at two different points yielding forward, central, and backward numerical differences respectively: \[ f_{\mathrm{FD}}' (x, h) := \frac{f(x+h) - f(x)}{h}, \quad f_{\mathrm{CD}}' (x, h) := \frac{f(x+h) - f(x-h)}{2h}, \quad f_{\mathrm{BD}}' (x, h) := \frac{f(x) - f(x-h)}{h} \]
These expressions involve only 2 terms, which limits their accuracy. Ignoring machine errors due to rounding, the approximation error of \(f_{\mathrm{FD}}' (x, h)\) and \(f_{\mathrm{BD}}' (x, h)\) due to Taylor series truncation is of the order \(O(h)\), and for \(f_{\mathrm{CD}}' (x, h)\), it is \(O(h^2)\); since in practice, \(h\) is very small, \(f_{\mathrm{CD}}' (x, h)\) is more accurate. The degree to which \(h\) is raised in the error term is reflected in the naming: \(f_{\mathrm{FD}}' (x, h)\) and \(f_{\mathrm{BD}}' (x, h)\) are first-order-accurate, while \(f_{\mathrm{CD}}' (x, h)\) is second-order-accurate. Higher accuracy can be achieved by computing more terms at extra points; the following approximation is fourth-order-accurate: \[ f_{\mathrm{CD},4}' (x, h) := \frac{f(x - 2h) -8 f(x-h) + 8 f(x+h) - f(x+2h)}{12h} \] Instead of \((x\pm h, x \pm 2h)\), any other convenient grid may be chosen. For the sake of brevity, we do not compare the accuracy of these approximations and refer to another vignette of this package, ‘Step-size-selection algorithm benchmark’.
Instead of choosing a large evaluation grid in advance, the researcher may instead compute the numerical derivative for several different step sizes, \(h_1 > h_2 > \ldots\), and combine multiple approximation to achieve higher accuracy. This is known as Richardson extrapolation: \[ f_{\mathrm{Rich},4}' (x, h_1, h_2) := \frac{(h_1 / h_2)^2 f_{\mathrm{CD}}' (x, h_2) - f_{\mathrm{CD}}' (x, h_1)}{(h_1 / h_2)^2 - 1} \]
If \(h_2 = h_1 / 2\), then, this extrapolation formula reduces to \(f_{\mathrm{CD},4}' (x, h_2)\). There is a general solution to obtaining high-order-accurate derivatives from function values evaluated at multiple points; see Fornberg (1988) for the solution and the algorithm and Kostyrka (2025) for derivation and stability analysis.
NB. For accuracy order \(a > 2\), more than 2 function evaluations are required, which is why the notion of ‘step size’ becomes unclear: it could signify the initial step size in the Richardson extrapolation before the first reduction, or the space between the elements of a symmetric grid \(\{\pm h, \pm 2h, \pm 3h, \ldots\}\) excluding zero, or the space between the elements of an equispaced symmetric grid \(\{\pm h, \pm 3h, \pm 5h, \ldots\}\), or some other measure describing the scaling of the deviations around the point of interest.
Therefore, to estimate a derivative numerically, the user should be able to supply the following tweaking parameters:
- The function \(f\) and evaluation point \(x\);
- The side to determine whether forward, central, or backward differences should be calculated;
- The step size \(h\);
- The desired accuracy order \(a\) for the truncation error to be \(O(h^a)\) (requiring at least \(2a\) evaluations). (Assuming, as described in the note above, that it means either the initial step size or the grid spacing based on the implementation.)
Vector cases and, therefore, Jacobians can be approximated by applying these formulæ to each coordinate of \(f\): \[ J_{f_4}(x) := \begin{bmatrix} \frac{\partial f^{(1)}_4}{\partial x^{(1)}} & \cdots & \frac{\partial f^{(1)}_4}{\partial x^{(n)}} \\ \vdots & \ddots & \vdots \\ \frac{\partial f^{(m)}_4}{\partial x^{(1)}} & \cdots & \frac{\partial f^{(m)}_4}{\partial x^{(n)}} \end{bmatrix}(x), \] where each partial derivative is estimated via finite differences.
Both numDeriv and pnd allow the user to
choose these parameters. The simplified comparison is given below.
| Argument | numDeriv::grad() syntax |
pnd::Grad() syntax |
|---|---|---|
| Function \(f\) | func |
FUN |
| Point \(x\) | x |
x |
| Side | side |
side |
| Step size | method.args = list(d = ..., eps = ...) |
h = ... |
| Accuracy order | method.args = list(r = ...) |
acc.order = ... |
1.1 Quick start: reproducing the numDeriv vignette
Examples of derivative and gradient calculation:
grad(sin, pi) # Old
#> [1] -1
Grad(sin, pi) # New
#> Estimated derivative: -1.0000
#> (default step size: 1.9e-05).Vector arguments are supported:
grad(sin, (0:10)*2*pi/10) # Old
#> [1] 1.000000 0.809017 0.309017 -0.309017 -0.809017 -1.000000 -0.809017
#> [8] -0.309017 0.309017 0.809017 1.000000
Grad(sin, (0:10)*2*pi/10) # New
#> Estimated gradient:
#> 1.0000 0.8090 0.3090 -0.3090 -0.8090 -1.0000 -0.8090 -0.3090 0.3090 0.8090 1.0000
#> (default step size: 1.5e-08, 5.2e-06, 7.6e-06, 1.1e-05, 1.5e-05, 1.9e-05, 2.3e-05, 2.7e-05, 3.0e-05, 3.4e-05, 3.8e-05).func0 <- function(x) sum(sin(x))
grad(func0, (0:10)*2*pi/10) # Old
#> [1] 1.000000 0.809017 0.309017 -0.309017 -0.809017 -1.000000 -0.809017
#> [8] -0.309017 0.309017 0.809017 1.000000
Grad(func0, (0:10)*2*pi/10) # New
#> Estimated gradient:
#> 1.0000 0.8090 0.3090 -0.3090 -0.8090 -1.0000 -0.8090 -0.3090 0.3090 0.8090 1.0000
#> (default step size: 1.5e-08, 5.2e-06, 7.6e-06, 1.1e-05, 1.5e-05, 1.9e-05, 2.3e-05, 2.7e-05, 3.0e-05, 3.4e-05, 3.8e-05).
x <- 2.04
numd1 <- grad(func1, x) # Old
numd2 <- Grad(func1, x) # New
numd3 <- Grad(func1, x, h = gradstep(func1, x)$par) # New auto-selection
exact <- 10*cos(10*x) + exp(-x)
c(Exact = exact, Old = numd1, New = numd2, NewAuto = numd2,
OldErr = (numd1-exact)/exact, NewErr = (numd2-exact)/exact,
NewAutoErr = (numd3-exact)/exact)
#> Exact Old New NewAuto OldErr NewErr NewAutoErr
#> 3.335e-01 3.335e-01 3.335e-01 3.335e-01 8.469e-12 -1.530e-09 2.977e-09Here, the relative error of the new approach with default settings is
worse than that of the old one. Nevertheless, the new implementation is
much faster and attains satisfactory accuracy with only 2 function
calls, whilst the numDeriv counterpart calls the function 9
times. If an appropriate step size is chosen (e.g. with the showcased
auto-selection procedure), then, the reliability of the result increases
while the derivative calculation still requires only 2 evaluations.
x <- c(1:10)
numd1 <- grad(func1, x) # Old
numd2 <- Grad(func1, x) # New
numd3 <- Grad(func1, x, h = sapply(x, function(y) gradstep(func1, y)$par))
exact <- 10*cos(10*x) + exp(-x)
cbind(Exact = exact, Old = numd1, New = numd2, NewAuto = numd2,
OldErr = (numd1-exact)/exact, NewErr = (numd2-exact)/exact,
NewAutoErr = (numd3-exact)/exact)
#> Exact Old New NewAuto OldErr NewErr NewAutoErr
#> [1,] -8.023 -8.023 -8.023 -8.023 6.118e-12 -6.443e-10 1.965e-09
#> [2,] 4.216 4.216 4.216 4.216 6.271e-12 -2.370e-09 2.843e-09
#> [3,] 1.592 1.592 1.592 1.592 -3.507e-12 -5.323e-09 2.863e-09
#> [4,] -6.651 -6.651 -6.651 -6.651 6.698e-12 -9.812e-09 1.502e-08
#> [5,] 9.656 9.656 9.656 9.656 -7.892e-13 -1.527e-08 1.344e-08
#> [6,] -9.522 -9.522 -9.522 -9.522 -4.610e-12 -2.200e-08 -1.448e-08
#> [7,] 6.334 6.334 6.334 6.334 1.243e-11 -2.993e-08 9.228e-09
#> [8,] -1.104 -1.104 -1.104 -1.104 6.093e-12 -3.913e-08 1.172e-08
#> [9,] -4.481 -4.481 -4.481 -4.481 1.524e-12 -4.952e-08 7.520e-09
#> [10,] 8.623 8.623 8.623 8.623 -8.120e-13 -6.111e-08 -2.153e-08Examples of Jacobian calculation:
func2 <- function(x) c(sin(x), cos(x))
x <- (0:1)*2*pi
jacobian(func2, x)
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
#> [3,] 0 0
#> [4,] 0 0
Jacobian(func2, x)
#> Estimated Jacobian:
#> 1.0000 1.0000
#> 0 0
#> (default step size range: 1.5e-08...3.8e-05.)Examples of Hessian calculation:
x <- 0.25 * pi
hessian(sin, x)
#> [,1]
#> [1,] -0.7071068
fun1e <- function(x) sum(exp(2*x))
x <- c(1, 3, 5)
hessian(fun1e, x, method.args=list(d=0.01))
#> [,1] [,2] [,3]
#> [1,] 2.955622e+01 2.583236e-15 1.620102e-11
#> [2,] 2.583236e-15 1.613715e+03 5.400691e-12
#> [3,] 1.620102e-11 5.400691e-12 8.810586e+041.2 Vectorisation pitfalls
Vectorisation does not occur naturally in computer science; it is rather a feature of a specially designed system than a natural property of hardware or software. From the theoretical point of view, computers are imperfect, restricted versions of Turing machine, analogous to deterministic finite automata or linear bounded automata, capable of reading and writing on a ‘tape’ of practically limited size provided a sequence of instructions. At the hardware level, computers are capable of executing a set of processor instructions, and these instructions are typically supplied sequentially. In fact, vectorisation at the hardware level is a relatively rare phenomenon because without a proper set-up, it may lead to race conditions and non-deterministic behaviour. In certain architectures or extensions (MMX, SSE, AVX), the data for calculations can be set up in such a manner that one single instruction processes multiple inputs, which can be achieved by packing the values into one large register (e.g. four 32-bit numbers into one 128-bit SIMD register). However, it is not easy to organise control flow in such a manner, which is why most computer operations are carried out sequentially.
The strong point of the R language is its
vectorisation. One may argue that this feature leads to its inefficiency
compared to, e.g., C++. Indeed, in C++, the overhead in loops is
typically irreducible and boils down to such hard-to-control factors as
jump instructions and data alignment to memory pages, whereas in R, the
overhead appear at each loop iteration due to environment manipulation,
extra assignment, name look-up, and occasional class-specific method
execution. On the other hand, not writing loops is useful when it would
take more human time to allocate memory and run a loop than computer
time to repeat redundant high-level operations under the hood: what
matters is the total time saved, machine and human combined. Finally, if
there is no explicit loop in a function because the latter relies
internally on apply() or calls C/C++ loops on high-level or
computationally heavy objects, then, the overhead due to the high-level
nature of R does not create substantial time losses compared to the
inner computations, and a vectorised function can be as slow as its
loop-based sequential version. This can be seen in the definition of
do_vapply in apply.c: the user-supplied
function is called in a loop; the implementation of
SEXP lapply also contains a C loop. Therefore, in an
efficient implementation of a vectorised function, the overhead is
minimised if most of the loops in which the input elements are processed
are low-level C/C++ loops with as few exchanges and high-level syntactic
embellishments as possible.
In R, it is harder to find something that is
not a vector. A scalar is just a special case of a vector –
with length 1: identical(c(1), 1) is TRUE.
Therefore, if a function is implemented in such a manner that it can
produce vector outputs for vector inputs en masse, it is
considered to be vectorised – with a caveat: dimensionality matters.
Consider the following classification:
| Scalar output | Vector output | |
|---|---|---|
| Scalar input | \(f_1\colon \mathbb{R} \mapsto \mathbb{R}\) | \(f_2\colon \mathbb{R} \mapsto \mathbb{R}^m\) |
| Vector input | \(f_3\colon \mathbb{R}^n \mapsto \mathbb{R}\) | \(f_4\colon \mathbb{R}^n \mapsto \mathbb{R}^m\) |
To compute a numerical derivative via finite differences, the function must be evaluated at least twice. This means that, depending on the parallel capabilities of a function, an efficient programmer should implement numerical differentiation in such a manner that functions that are vectorised at the low level take as many arguments in one call as possible. Many calls of the same function with shorter arguments are slower than one call with a long argument:
system.time(replicate(1e5, sin(rnorm(10))))
#> user system elapsed
#> 0.271 0.008 0.278
system.time(sin(rnorm(1e6)))
#> user system elapsed
#> 0.073 0.000 0.074The same phenomenon is observed in the real world with storage media: copying 1000 files of size 1 MB each onto a flash drive is slower than copying one 1-GB file obtained by concatenating said 1000 files.
The case of non-vectorised scalar output for scalar input, \(f_1\), can be vectorised trivially by
putting the function inside a loop or an apply-like
operation. Then, the preparation task for approximating \(f'_1\) is creating an array of values
at which the function will be evaluated in a loop or in a sequence of
parallel jobs (if \(f_1\) is
computationally costly and the parallelisation overhead is small
compared to the evaluation time of \(f_1\)). This is the case if \(f_1(x_0)\) has length 1 for scalar \(x_0\) and throws and error for vector \(x_0\).
If the function \(f_2\) returns a vector of length \(m\) for scalar input, then, its coordinate-wise derivatives make up a Jacobian with dimensions \(m\times 1\). Such a function should return an error for vector input, but a list of scalar inputs should produce a list of vector outputs in a loop / apply-like call that can be combined to obtain the Jacobian approximation. The user should be aware of this function behaviour and specifically request a Jacobian.
The case of \(f_3\) is simple: if the output has length 1 for vector input of length \(n\), then it implies that \(f_3\) carries out a dimensionality-reducing operation on the input elements. This implies that the numerical derivatives with respect to each coordinate form the gradient vector, and the user should call the function that computes the gradient approximation. The numerical derivative with respect to the first coordinate of \(x_0\) is the weighted sum of \(f(x_0^{(1)} + h, x_0^{(2)}, \ldots)\), \(f(x_0^{(1)} - h, x_0^{(2)}, \ldots)\), and potentially other values obtained at points where only the first coordinate of \(x_0\) is altered. This is repeated for each input coordinate: e.g., for central differences, \(2n\) evaluations are expected. The input is therefore a list of \(2n\) argument vectors of length \(n\), and potentially more elements for higher accuracy or higher derivation order.
The case with \(f_4\) is the most
complicated one, and it is the reason behind the unexpected behaviour of
some of the functions in numDeriv. There are two sources of
danger in this implementation.
1. If \(n = m\), i.e. the user supplies a vector of length \(n\) and the function returns a vector of length \(n\), it is either due to the true vectorisation of \(f_4\) or due to the fact that the output always has length \(n\) regardless of the input, and such a coincidence happened that the input had length \(n\). It is impossible to claim that \(f_4\colon \mathbb{R}^n \mapsto \mathbb{R}^n\) is a vectorised \(\vec f_1\colon \mathbb{R} \mapsto \mathbb{R}\) solely from the fact that the output has length \(n\) for some input of length \(n\); multiple lengths should be checked to confirm this.
Example 1. Consider a function that computes the
quartiles of its input series. Its output always has length 3.
numDeriv::grad will believe that this function is a
vectorised scalar-valued one (i.e. \(\vec
f_1\)) if length(x) == 3:
f <- function(x) quantile(x, 1:3/4)
grad(f, x = 1:2)
#> Error in grad.default(f, x = 1:2): grad assumes a scalar valued function.
grad(f, x = 1:4)
#> Error in grad.default(f, x = 1:4): grad assumes a scalar valued function.
grad(f, x = 1:3)
#> [1] 1.5000000 1.0000000 0.8333333It is therefore necessary to explicitly call its Jacobian:
jacobian(f, x = 1:3)
#> [,1] [,2] [,3]
#> [1,] 0.5 0.5 0.0
#> [2,] 0.0 1.0 0.0
#> [3,] 0.0 0.5 0.5
jacobian(f, x = 1:4)
#> [,1] [,2] [,3] [,4]
#> [1,] 0.25 0.75 0.00 0.00
#> [2,] 0.00 0.50 0.50 0.00
#> [3,] 0.00 0.00 0.75 0.25This occurrence is more likely if \(n\) and \(m\) are small: for long inputs, a long output vector of identical length implies almost certainly that \(f_4\) is a vectorised \(\vec f_1\colon \mathbb{R} \mapsto \mathbb{R}\), enabling efficient evaluation grid preparation to keep the overhead due to loops as low as possible. \(\blacksquare\)
2. If \(n\ne m\) but the function is internally parallelised coordinate-wise, then, dimension checks based on the observed input and output lengths may be flawed and result in unexpected behaviour. Technically, computing \(f_4\colon \mathbb{R}^n \mapsto \mathbb{R}^m\) can be seen as repeated application of \(f_2\colon \mathbb{R} \mapsto \mathbb{R}^m\) to a list of \(n\) points, and its parallelisation can be envisioned as such. However, if \(f_4\) is implemented as stacked \(f_3\colon \mathbb{R}^n \mapsto \mathbb{R}\), the output should be trated differently.
The Jacobian matrix can be created by column, stacking evaluations of \(\frac{\partial}{\partial x} f^{(i)}_4(x)\) at each point \(x_1, \ldots, x_n\), or by row, binding the transposed gradients \(\nabla^\intercal f^{(i)}_4(x)\) of each coordinate of \(f\). In R, matrices can be created by column (default) or by row, which determines the order in which the vector of length \(nm\) is wrapped, but the exact method has to be decided by the user or, in case of a package, inferred from the function behaviour.
Example 2. Consider \(f(x)
:= \begin{pmatrix} \sin x \\ \cos x \end{pmatrix}\). Note that
the functions sin() and cos() are vectorised
and can handle inputs of arbitrary lengths.
f <- function(x) c(sin(x), cos(x))
f(1)
#> [1] 0.8414710 0.5403023
jacobian(f, 1:2)
#> [,1] [,2]
#> [1,] 0.5403023 0.0000000
#> [2,] 0.0000000 -0.4161468
#> [3,] -0.8414710 0.0000000
#> [4,] 0.0000000 -0.9092974
jacobian(f, 1:3)
#> [,1] [,2] [,3]
#> [1,] 0.5403023 0.0000000 0.0000000
#> [2,] 0.0000000 -0.4161468 0.0000000
#> [3,] 0.0000000 0.0000000 -0.9899925
#> [4,] -0.8414710 0.0000000 0.0000000
#> [5,] 0.0000000 -0.9092974 0.0000000
#> [6,] 0.0000000 0.0000000 -0.1411200The Jacobian of \(f\) is \((\cos x, -\sin x)\); therefore, the
expected output for x0 = 1:2 is \(\begin{pmatrix} \cos 1 & \cos 2 \\ -\sin 1
& -\sin 2 \end{pmatrix}\). However, the output is a matrix
made of stacked diagonal matrices, and the result is wrong even when
\(n\ne m\)!
The root cause of this error is the manner in which the vector-valued
function is defined: the length of the object
c(sin(x), cos(x)) depends on the length of x:
if length(x) == n, then,
length(c(sin(x), cos(x))) == 2*n, and R
sees the output as a vector, not a matrix! Hence,
jacobian(f, x0) computes f(x0) to check the
function dimension, infers from the output
f(c(sin(1:2), cos(1:2))) that \(m
= 4\), \(n = 2\), which is
wrong, and pre-allocates a \(4\times
2\) matrix of zeros for the results. \(\blacksquare\)
Therefore, without any a priori information about the function inputs and outputs, the function that computes matrices of derivatives should follow the following logic for determining dimensions and handling potential errors:
- If the input argument
xis a scalar, prepare an evaluation grid in the form of a list, evaluatefon that grid, and combine the result by row. For safety, ignore the fact thatfmay be vectorised (likesin(1:100)) and may handle vector inputs in some cases.- If the resulting matrix has 1 column, then, the function is either
\(f_1\colon \mathbb{R} \mapsto
\mathbb{R}\) or \(f_3\colon
\mathbb{R}^n \mapsto \mathbb{R}\) (since the output has length
1). Proceed with computing the derivative approximation via weighted
sums and return the derivative estimate. If a Jacobian was initially
requested, throw a warning that the output of
fis a scalar and gradient facilities should be used. - If the resulting matrix has more than 1 column, then, the function
is either \(f_2\colon \mathbb{R} \mapsto
\mathbb{R}^m\) or \(f_4\colon
\mathbb{R}^n \mapsto \mathbb{R}^m\). Proceed with computing the
Jacobian approximation via weighted sums and return the row matrix of
the Jacobian estimate. If a gradient was initially requested, throw a
warning that the output of
fis not a scalar and Jacobian facilities should be used, or return an error that this case should be handled by the Jacobian routine directly.
- If the resulting matrix has 1 column, then, the function is either
\(f_1\colon \mathbb{R} \mapsto
\mathbb{R}\) or \(f_3\colon
\mathbb{R}^n \mapsto \mathbb{R}\) (since the output has length
1). Proceed with computing the derivative approximation via weighted
sums and return the derivative estimate. If a Jacobian was initially
requested, throw a warning that the output of
- If the input argument
xis a vector, first, it must be determined iffis vectorised \(\vec f_1\colon \mathbb{R} \mapsto \mathbb{R}\) or not. Callf(x)and examine the output.- If the output has the same length as the input, it is highly likely
that
fis vectorised – for a counterexample where it is not, see Example 1 above showing that if the output is a vector, an input of identical length may create the impression of vectorisation by pure chance. This case can be checked by callingf(x[1])and checking if the output has length 1. If not, a warning that the user should request a Jacobian for a vector-valued function should be issued. If this check is passed, proceed by assuming vectorisation and prepare the evaluation grid as a long vector to reduce overhead. - If the output has length 1, then, the function is certainly \(f_3\colon \mathbb{R}^n \mapsto
\mathbb{R}\), and
fshould be applied to a list of input vectors for gradient approximation. If the user requested a Jacobian, a warning should be shown. - If the output length is neither
length(x)nor1, thenfdefinitely belongs to the class \(f_4\colon \mathbb{R}^n \mapsto \mathbb{R}^m\). In this case, similarly to the the case with scalar output, a list of vectors should be prepared as the evaluation grid for Jacobian approximation. If the user requested a gradient, throw a warning that the result is a Jacobian matrix with non-null dimensions.
- If the output has the same length as the input, it is highly likely
that
The user should know whether their function is scalar- or vector-valued and call the gradient or Jacobian routine accordingly. Nevertheless, a check should be carried out for dimensionality and vectorisation.
Gradient vectorisation detection is implemented as
follows. Firstly, compute f(x).
- If
length(x) == 1andlength(f(x)) == 1, create the evaluation grid as a list and useapply-like calls because there is no knowledge whetherfis \(f_1\) or \(\vec f_1\). - If
length(x) > 1andlength(f(x)) == 1, create the evaluation grid as a list and useapply-like calls. This is the \(f_3\) case. - If
length(x) > 1andlength(f(x)) > 1, check iflength(f(x)) == length(x), create the evaluation grid as a vector and callfon it directly (the \(\vec f_1\) case). Otherwise, return an error that a Jacobian should be computed instead because the input vector length is not equal to the output vector length. - If
length(x) == 1andlength(f(x)) > 1, return an error indicating that a Jacobian should be computed instead.- Checking if
length(f(x[1])) == 1to avoid the error from Example 2 above is unreliable because certain functions (especially statistical ones) are not defined for scalars and might return an error. Example: a function that returns the deviations from the trimmed mean where the highest and the lowest values are forcibly deleted could fail if the number of observations is less than 3. Similarly, if the input length determined the number of observations for a model, having a single observation may result in non-invertibility of the model matrix; this error would not occur for long vector inputs.
- Checking if
Jacobians should be checked similarly: compute
f(x) and check the following.
- If
length(f(x)) == 1, return an error indicating that a gradient should be computed instead. - If
length(x) > 1andlength(f(x)) > 1, create the evaluation grid as a list and callfon it directly. This is the \(\vec f_4\) case. Even if the components of \(f\) are vectorised, applyingfto a list of points is the only safe option.
