10
$\begingroup$

What is the best method for 1D numeric differentiation? Something as glorious as Gaussian quadrature for numeric integration.

Maybe differential quadrature is such a method? What is its accuracy?

I'm well aware that it is really easy to have symbolic differentiation in the program (automatic differentiation or truly symbolic algorithm). However to use such methods it is necessary to rewrite all functions to be differentiated. Thus one can't differentiate functions imported from libraries.

I need differentiation almost with the machine precision.

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ J.M. has already supplied some excellent ways of performing numerical differentiation, including methods like Cauchy's formula. There is a simple method that I use that gets me first derivatives at near machine precision levels. It is the complex step derivative: $f'(x) = \mathrm{Im}(\frac{x + ih}{h})$, where $h$ can be chosen to be the machine epsilon (see citeseerx.ist.psu.edu/viewdoc/… for guidance). It is trivial to implement in any language with a complex number datatype (e.g. Fortran). $\endgroup$
    – Gilead
    Commented May 8, 2011 at 16:23
  • $\begingroup$ Ok, I see that J.M. has linked to this method (due to Squire and Trapp) in one of the references. $\endgroup$
    – Gilead
    Commented May 8, 2011 at 16:24
  • 1
    $\begingroup$ @Gilead: The only reason why I linked to it in the comments instead of my answer is that even though I know it's good, I haven't extensively experimented with it. $\endgroup$
    – J. M. isn't a mathematician
    Commented May 8, 2011 at 17:07
  • 1
    $\begingroup$ BTW... I trust that you only want first derivatives? The higher the derivative order, the less likely any of the proposed methods become successful... $\endgroup$
    – J. M. isn't a mathematician
    Commented May 8, 2011 at 18:03
  • 9
    $\begingroup$ It's worth noting that the complex step derivative is just a thinly disguised version of automatic differentiation. It works because real analytic functions have their derivatives baked into their complex implementations. $\endgroup$
    – Geoffrey Irving
    Commented Mar 10, 2017 at 3:47

5 Answers 5

Reset to default
18
$\begingroup$

If your function is badly behaved (e.g. noisy, very oscillatory), no method will perform properly (differentiation is numerically very unstable). That being said, for "nice functions", I have good experience with polynomial (Richardson) extrapolation methods. This paper and this paper give hints on how you might write your own implementation. I will note that this is the method implemented in the NAG numerical libraries (with of course a few wrinkles of their own).

There are two possible alternatives if for some reason you don't want to use Richardsonian methods. One is to use Cauchy's differentiation formula:

$$f^\prime(x)=\frac1{2\pi i}\oint_\gamma \frac{f(t)}{(t-x)^2}\mathrm dt$$

where it is up to you to choose a suitable counterclockwise contour $\gamma$ (a circle is customary); the other is to use the "Lanczos derivative":

$$f^\prime(x)=\lim_{h\to 0}\frac{3}{2h^3}\int_{-h}^h t\;f(x+t)\mathrm dt$$

where you either will have to experiment with an appropriate step size $h$, or use some extrapolative procedure.

You will have to experiment with your computing environment to choose.

Improve this answer
$\endgroup$
9
  • 2
    $\begingroup$ This algorithm from the ACM might also be of interest: dx.doi.org/10.1145/362759.362820 (@Netlib: netlib.org/toms/413). See also dx.doi.org/10.1137/0704019 , ams.org/journals/mcom/1968-22-102/S0025-5718-1968-0230468-5/… , dx.doi.org/10.1137/S003614459631241X , and dx.doi.org/10.1145/838250.838251 . You've quite a lot to choose from, if you are able to evaluate your function to differentiate at complex arguments (as with using the Cauchy differentiation formula). $\endgroup$
    – J. M. isn't a mathematician
    Commented May 8, 2011 at 16:04
  • $\begingroup$ Thank you very much for the respond, now I'm studying the techniques connected to Lanczos derivative, at least at the moment they look promising for me. I will report later on the results. As for using complex numbers I don't actually see the point in it. I can't use them to find derivatives of already implemented functions of type $\text{Double} \to \text{Double}$. If I was to reimplement these functions I would rather use dual numbers to have exact (i.e. exact as symbolical ones) derivatives. $\endgroup$
    – Yrogirg
    Commented May 8, 2011 at 17:41
  • 1
    $\begingroup$ "As for using complex numbers I don't actually see the point in it." - the differentiation problem tends to be a bit more stable if you do complex evaluation versus constraining yourself to evaluating only at real values, but that's why I gave the Richardson and Lanczos methods since I do know that sometimes you can only evaluate at reals... $\endgroup$
    – J. M. isn't a mathematician
    Commented May 8, 2011 at 17:57
  • 1
    $\begingroup$ At the end I have abandoned all the numeric techniques and decided to stick to symbolic differentiation. I had to rewrite the program but now at least I have exact differentiation. May be later I will have to return to the issue. $\endgroup$
    – Yrogirg
    Commented May 10, 2011 at 15:18
  • 5
    $\begingroup$ I don't understand how the Lanczos derivative is any better than the finite-difference approximation. In order to evaluate the integral you need to perform a subtraction, which seems to be exactly equivalent to the finite-difference formula in terms of accuracy/stability... am I missing something? How are you supposed to do that integration? $\endgroup$
    – user541686
    Commented Mar 1, 2015 at 13:02
10
$\begingroup$

Use auto-differentiation. Automatic differentiation is faster than other forms of differentiation and gets errors at machine epsilon. However, it's much harder to implement, so you might need to use a package for it. That said, since it's so useful there are plenty of packages, one I use very often is ForwardDiff.jl.

It's a shame that less academics don't use automatic differentiation because it really is what should be the workhorse, and it's just not taught in graduate school for some reason (though those who go into machine learning know it as "backwards propagation").

Edit

One can autodifferentiate functions from libraries in Julia since the libraries themselves are written in Julia and are "duck-typed", and so ForwardDiff.jl can auto-differentiate library functions (including ones which use loops, conditionals, etc.) due to its type system. The only problem it runs into is when there are C/Fortran libraries called, so really you just have to avoid library functions which internally call BLAS (so no matrix multiplications and you're good).

Edit 2

ForwardDiff will now use generic linear algebra fallback function to compute the derivatives, so those will be fine as well. The only Base library functions I can think of where automatic differentiation will fail now (in Julia) is when trying to differentiate a function which calculates an FFT, since that uses FFTW and doesn't have a Julia fallback (for now).

Improve this answer
$\endgroup$
2
  • 5
    $\begingroup$ I agree that automatic differentiation is great and underutilized, but the OP explains why it doesn't apply in the case at hand. $\endgroup$
    – Noah Stein
    Commented Nov 15, 2016 at 1:13
  • 1
    $\begingroup$ I forget to add the whole reason why I posted this. It does apply to this situation if you're using the right tools. I made sure to note that. $\endgroup$
    – Chris Rackauckas
    Commented Nov 15, 2016 at 1:22
5
$\begingroup$

Complex step differentiation (CSD) is well known as an efficient numerical differentiation method: $$f^\prime(x)=\mathrm{Im}\frac{f(x+\mathrm{i}h)}{h}+O(h^2),\quad\mathrm{i}:=\sqrt{-1}.$$ This method requires the function to be analytic (differentiable as a complex function). We cannot apply this method to higher order derivatives, but since this method does not use differences it is useful to avoid subtractive cancellation (loss of significance caused when we try to subtract similar numbers). For the derivation and implementations, check the article on SIAM News and the post on Mathworks Blog.

A noteworthy alternative to CSD is the approach by dual numbers. Dual numbers are defined by: $$x=a+b\epsilon,\quad a,b\in\mathbb{R},\quad\epsilon^2=0.$$ By using Taylor expansion and dual numbers, we obtain $$f(a+b\epsilon)=f(a)+bf^\prime(a)\epsilon\quad(\because\epsilon^2=0).$$ Hence we can obtain the derivative by dual numbers.

Improve this answer
$\endgroup$
5
$\begingroup$

In addition to J.M. answer: this paper https://www.sciencedirect.com/science/article/pii/S0021904512000123 (Differentiation by integration using orthogonal polynomials, a survey, by E. Diekema and T.H. Koornwinder, doi: 10.1016/j.jat.2012.01.003) provides further insight into Lánczos's derivative and its generalizations for n-th order derivatives. By the way, it says that the Lánczos's (1956) derivative formula goes back to Cioranescu (1938) and Haslam-Jones (1953).

Improve this answer
$\endgroup$
2
$\begingroup$

Warning: self-promotion

I presume the "numerical" aspect of numerical differentiation is very well studied, what is little bit forgotten is "statistical" aspect. To increase the precision of numerical differentiation do the following:

1) Chose your favorite high-precision "standard" method based on some step size H.

2) Compute the value of the derivative with the method chosen in 1) many times with different but reasonable step sizes h. Each time you may pick h as a random number from the interval (0.5*H/10, 1.5*H/10) where H is an appropriate step size for the method you use.

3) Average the results.

Your result may gain several orders of magnitude in the absolute error wrt. the non-averaged result. It is time-consuming method, but we discuss here precision and not time consumption. I franky believe it is the most precise method of numerical differentation at the market today (yet systematically refused by journals)

https://arxiv.org/abs/1706.10219

Improve this answer
$\endgroup$
8
  • $\begingroup$ The answer would be vastly improved by giving a summary of the contents of the link. $\endgroup$
    – Tommi
    Commented May 18, 2018 at 9:35
  • 2
    $\begingroup$ "we discuss here precision and not time consumption" This makes no sense - the strength of a method is basically how much time is needed to achieve a given precision (except when memory complexity is the limiting factor, but I do not see a reason to expect this to be the case here). $\endgroup$
    – Benoît Kloeckner
    Commented May 20, 2018 at 21:00
  • 1
    $\begingroup$ You are wrong to say "strength" of the method. The question clearly says "precision" and this is discuddes here. As mentioned: one downloads numerical libraries he did not program and one wants the most precision when differentiating: it is a well defined problem. $\endgroup$
    – F. Jatpil
    Commented May 20, 2018 at 21:34
  • 1
    $\begingroup$ From my point of view I can see why it is "systematically refused by journals". There is no guarantee that 'averaging out' rounding errors will yield better results, neither does this improve convergence order or similar. I can get reliably get higher accuracy by increasing my precision, with comparable increase of computational effort. $\endgroup$
    – oli
    Commented May 27, 2020 at 20:14
  • 1
    $\begingroup$ @F.Jatpil: In your linked paper the tables imply that if h is close to properly chosen for the specific method, then you get 1 or 2 digits at most of increased accuracy. In the case of central finite differences I expect this improvement to vanish if h is exactly properly chosen and rounding errors are avoided for (x+h) and (x-h). If you wanna argue that we can perturb floating arithmetic calculations several times, average the result, and maybe get one digit or two in accuracy, fine, that might often work. I don't see a real value in that. $\endgroup$
    – oli
    Commented May 28, 2020 at 13:56

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .