Computing multiple numerical derivatives at once

Question

Lets say I have a function $f(X) = f(x_1,...,x_N)$ to be integrated. But unlike time discrete methods, my integrator uses quantisation to advance time, that is if $|x - q| > dQ$, with $q$ being the quantised version of $x$ and $dQ$ being the quantum.

To achieve this, the function needs to send $f(X)$, $f'(X)$, and $f''(x)$ to the integrator, which determines the time step, $t = \sqrt[3]{( dQ \frac{3}{f''(x)} )}$, from the 2nd derivative (of the derivative function $f$) and after $t$ units time sends back the updated (integrated) value ($[v,mv,pv]$ ~ value, 1st & 2nd quantized derivative values of $f$).

In the paper, $f'$ is approximated by calculating

$c_j = [ f(X) - f(X^*)]/ (x_j - x^*_j)$ for each updated value $x^*_j$

and then

$f' = \sum( c_j ) \cdot mv$,

But in the reference implementation, they update all values, and then simply calculate

$f' = f(X) - f(X^*) / 10^{-16} $(and so on for the 2nd derivative)

I'm unsure on which way would be preferable, the first approach using $2m$ evaluations, the second one just calls $f$ twice.

Which way would be better for non-linear function approximation?, given that it needs to be accurate only within the quantum bounds (which are generally small,. e.g $0.0001 \cdot |f(X)|$ )

Hint: this is a discrete event method, large jumps might be present.

Answer 1

It's not entirely clear to me what you're asking. Let me see answer what I think you're asking, and you can tell me if I'm on the right track.

Suppose you've got a function $f: \mathbb{R}^{m} \rightarrow \mathbb{R}$, which is to say that it takes $m$ variables $x_{1}, \ldots, x_{m}$ as inputs, and outputs a scalar. Also, suppose you want to evaluate derivatives of $f$ at some point $\bar{\mathbf{x}} = (\bar{x}_{1}, \ldots, \bar{x}_{m})$.

If you want to evaluate $\frac{\partial{f}}{\partial{x_{i}}}(\bar{\mathbf{x}})$, then you could approximate it by

\begin{align} \frac{f(\bar{\mathbf{x}} + h\mathbf{e}_{i}) - f(\bar{\mathbf{x}})}{h}, \end{align}

where $\mathbf{e}_{i}$ is a unit vector of all zeros, except for a 1 in the $i$th position.

More generally, if you want to evaluate $\nabla{f}(\bar{\mathbf{x}}) \cdot \mathbf{v}$, you could approximate it by

\begin{align} \frac{f(\bar{\mathbf{x}} + h\mathbf{v}) - f(\bar{\mathbf{x}})}{h}. \end{align}

In both cases, you will want to use a small $h$; how small depends on the application. Good guidelines for choosing $h$ can be found in the review paper Jacobian-free Newton-Krylov methods: a survey of approaches and applications.

The finite difference approximation does degrade somewhat the accuracy of the derivative, so if your calculations require precise derivative information, you are better off trying other methods (depending on application, these could be forward or adjoint sensitivity methods; automatic differentiation; complex number finite differences; using some combination of pencil, paper, and your favorite computer algebra system) that get you exact (or close to exact) derivative information.

Applying this process twice to get a Hessian only further degrades the quality of the approximation, and requires $m^{2}$ operations, and is generally regarded as too expensive and too inaccurate to be worth doing. Applications that do need Hessian information can sometimes get by with approximation Hessian information obtained from an initial guess plus rank-1 updates to this approximate Hessian using gradient information. One example would be BFGS methods.