why is complex step differentation better than central differences

I came across this post at Cleve's Corner, where he shows that complex step differentiation is more accurate than central differences. The error in both methods is $$O(h^2)$$. So why exactly does complex step differentiation yield much better results? Is is because the constants hidden in $$O(h^2)$$ are much better for complex step differentiation?

• This is due to cancelation of significant digits in the difference of the numerator, which does not occur in the complex step differentiation. See e.g. Nick Higham's talk here.
– Bort
May 14 at 6:14
• Wow, that's clever. Has this technique been used for simulations involving solving time-evolution PDEs by Newton-Krylov, e.g., CFD? May 14 at 7:37
• I think I see the source of my confusion. The methods have the same truncation error but not the same roundoff error May 14 at 8:34
• You could answer your own question. That way somebody will find it later. May 14 at 16:43
• @MaximUmansky: there is a paper from the 2000's where they apply the CSD to solve a PDE. However, due to the fact nobody uses it nowadays, it doesn't seem to be advantageous for PDEs. In this answer, I've used the CSD to calculate Finite-difference coefficients on arbitrary grids, which I personally find quite nice. May 15 at 0:04

Indeed, the roundoff error is much lower with the complex step approach, because all the derivative terms are handled in the imaginary part, without direct interaction with the real part, as opposed to a finite difference calculation.

The advantage of complex step is that you can take a very low perturbation size $$h$$ (I usually use $$10^{-50}$$), whereas you would only go down to around $$h=\sqrt{\epsilon_{machine}} \approx 10^{-8}$$ for a classical finite difference formula to avoid excessive roundoff errors (see the convergence plots in your linked Matlab post for instance).

Indeed, let's consider the Taylor series of $$f$$ with a complex perturbation: $$f(x+ih) = f(x) + ihf'(x) - h^2 f''(x) - ih^3f'''(x) + O(h^4) + i O(h^5)$$

Taking the imaginary part and dividing by $$h$$, you get: $$f'(x) - h^2 f'''(x) + O(h^4)$$ By taking a very small $$h$$, the second and third term completely disappear due to roundoff errors in the imaginary part, thus leaving you with a machine-precision accurate evaluation of $$f'(x)$$, which you could most likely never attain with a classical finite difference scheme (or maybe with a very high-order scheme). An alternative would be to use automatic differentiation.

Another aspect is that the complex step is usually easy to implement (if the language allows for handling complex variables easily). For instance for a scalar function $$f$$, you can construct its first derivative in Python as a one-liner:

fprime = lambda x: np.imag( f(x + 1e-50*1i) ) * 1e50

• This is very true and useful assuming that the function is analytic at the point of interest. However, the first order FD approximation would work as long as the function is differentiable, so it is robust in a different sense. I don't have a pathologic enough function off the top of my head for which CD would fail, but FD would succeed; but there are indeed smooth but non-analytic functions of single variable (see the Wikipedia page en.wikipedia.org/wiki/Non-analytic_smooth_function). May 21 at 8:32
• Indeed, there might be some cases where CS does not work. Yet, I am still to find one in my use cases ! I would be curious to know if you find an applied case that suffers from this. May 21 at 13:47
• I suspect I will not find an example in practice; but that will be mostly because I usually don't use CS and whenever I use it I usually have the assumption that the function is analytic, so there is no reason to double-check. But I found an exercise from a complex analysis textbook: $f(z)$ is piecewise defined with $f(z) = z^5/\|z^4\|$ if $z\neq 0$ and $f(z) = 0$ otherwise. At $z=0$, the partial derivative in the direction of the real axis is $1$, but the partial derivative in the direction of the imaginary axis is $i$. A solution like this may arise from solutions of fractional PDEs. May 21 at 15:24
• One further comment regarding automatic differentiation: the CSD actually is automatic differentiation in disguise, as the complex dimension is equivalent to the difference parameter in the AD forward mode (this works well in 1D, but in higher dimensions AD is the easier concept) May 22 at 9:06

Answering as suggested by @nicoguaro. The methods have the same truncation error but not the same roundoff error, that is why they have different numerical accuracy.