## Complex Magic, Part 1: Differentiation via Function Evaluation at a Single Point

At the heart of conventional methods of derivative approximation based on divided differences lies a sum of the form ${\sum a_i f(x_i)}$ in which the coefficients ${a_i}$ have different signs and vary (usually considerably) in magnitude. As a result, when these coefficients of different magnitude multiply function values, the absolute errors in the latter are amplified and the error in the outcome becomes quantitatively unpredictable [LynessMoler1967].

As a special case, consider the second-order central difference formula, one of the standard derivative approximation methods for a real valued function,

$\displaystyle f'(x_0) \approx (f(x_0+h) - f(x_0-h))/2h, \ \ \ \ \ (1)$

for which, as the name should suggest, the truncation error is ${O(h^2)}$. In theory, the truncation error could be made arbitrarily small by making the step size small enough. In practice, on a finite-precision machine with floating point arithmetic, this cannot be done due to the existence of another type of error, the rounding error. In the particular case of (1), for instance, too small a step size would lead to a catastrophic subtractive cancellation in the numerator, whose culprit is the rounding done when evaluating the exponential. Clearly, a compromise must be made between minimizing truncation and rounding errors. In the case of (1), analyzing the monotonic behavior reveals that when rounding becomes prominent, the approximate value given by (1) oscillates around the true value [SteplemanWinarsky1979]. Although this behavior can be used to select an optimal ${h}$, only marginal improvements can be made [SquireTrapp1998].

Using complex variables it is possible to derive derivative approximation formulas that do not involve subtractions and hence are free from catastrophic cancellations. Let us derive a complex-step numerical approximation formula for the first derivative with the same order of truncation error as (1). If ${f(z)}$ is an holomorphic function of a complex variable ${z}$ that is real valued on the real axis, then its Taylor expansion, with complex step ${ih}$, about the real point ${x_0}$ is

$\displaystyle f(x_0 + ih) = f(x_0) + ihf'(x_0) + (ih)^2 f''(x_0)/2! + (ih)^3 f'''(x_0)/3! + \cdots. \ \ \ \ \ (2)$

Taking imaginary parts of both sides of (2) yields

$\displaystyle f'(x_0) = \mathop{Im}[f(x_0 + ih)]/h + h^2 f'''(x_0)/3! - \cdots. \ \ \ \ \ (3)$

The first term on the RHS of (3) is an ${O(h^2)}$ estimate to the LHS which is the first derivative. As promised, ${\mathop{Im}[f(x_0 + ih)]/h}$ does not involve subtractions and, consequently, it will not suffer from the associated cancellation due to rounding.

This approach can be easily implemented to compute the jacobian matrix. (The implementation below is based on Yi Cao’s.)


function [J,z] = jacobianComplexStep(fun,x,varargin)

z = fun(x);
n = numel(x);
m = numel(z);
J = zeros(m,n);

h = n*eps;

optArgs = {h};
emptyArgs = cellfun(@isempty,varargin);
[optArgs{~emptyArgs}] = varargin{~emptyArgs};
[h] = optArgs{:};

for k = 1:n
x1 = x;
x1(k)  = x1(k) + h*1i;
J(:,k) = imag(fun(x1))/h;
end



Here is an equivalent implementation of the jacobian based on the second order central difference formula:


function [J,z] = jacobianCentralDiff(fun,x,varargin)

z = fun(x);
n = numel(x);
m = numel(z);
J = zeros(m,n);

h = n*eps;

optArgs = {h};
emptyArgs = cellfun(@isempty,varargin);
[optArgs{~emptyArgs}] = varargin{~emptyArgs};
[h] = optArgs{:};
h_2 = 2*h;

for k = 1:n
[x_left,x_right] = deal(x);
x_right(k)  = x(k) + h;
x_left( k)  = x(k) - h;
J(:,k)      = (fun(x_right) - fun(x_left))/h_2;
end



And here is a sample invocation of the functions jacobianComplexStep and jocobianCentralDiff:


f = @(x) x^(9/2);
x = 1.5;

le = -20;
ri = -2;
m  = ri - le + 1;
h  = logspace(ri,le,m).';

[dfdxComplexStep,dfdxCentralDiff] = deal(zeros(m,1));
for iRow=1:m
dfdxComplexStep(iRow) = jacobianComplexStep(f,x,h(iRow));
dfdxCentralDiff(iRow) = jacobianCentralDiff(f,x,h(iRow));
end

fprintf('%.16e %.16e %.16e\n',[h dfdxComplexStep dfdxCentralDiff]')



The tables below compare the two methods by listing their approximations to the first derivative of the function ${f(x) = x^{9/2}}$ at ${x_0 = 1.5}$ for different step sizes, using MATLAB R2009a on a machine running Windows XP with an Intel Core 2 Duo E8400 CPU.

 ${h}$ ${f'(1.5)}$, Complex Step ${1.0000000000000000\times 10^{-2}}$ ${1.8599607128036336\times 10^{1}}$ ${1.0000000000000000\times 10^{-3}}$ ${1.8600800678177631\times 10^{1}}$ ${1.0000000000000000\times 10^{-4}}$ ${1.8600812613698931\times 10^{1}}$ ${9.9999999999999991\times 10^{-6}}$ ${1.8600812733054148\times 10^{1}}$ ${9.9999999999999995\times 10^{-7}}$ ${1.8600812734247697\times 10^{1}}$ ${9.9999999999999995\times 10^{-8}}$ ${1.8600812734259637\times 10^{1}}$ ${1.0000000000000000\times 10^{-8}}$ ${1.8600812734259762\times 10^{1}}$ ${1.0000000000000001\times 10^{-9}}$ ${1.8600812734259762\times 10^{1}}$ ${1.0000000000000000\times 10^{-10}}$ ${1.8600812734259762\times 10^{1}}$ ${1.0000000000000001\times 10^{-11}}$ ${1.8600812734259758\times 10^{1}}$ ${9.9999999999999998\times 10^{-13}}$ ${1.8600812734259762\times 10^{1}}$ ${1.0000000000000000\times 10^{-13}}$ ${1.8600812734259765\times 10^{1}}$ ${1.0000000000000000\times 10^{-14}}$ ${1.8600812734259762\times 10^{1}}$ ${1.0000000000000001\times 10^{-15}}$ ${1.8600812734259762\times 10^{1}}$ ${9.9999999999999998\times 10^{-17}}$ ${1.8600812734259762\times 10^{1}}$ ${9.9999999999999992\times 10^{-18}}$ ${1.8600812734259758\times 10^{1}}$ ${1.0000000000000001\times 10^{-18}}$ ${1.8600812734259762\times 10^{1}}$ ${9.9999999999999998\times 10^{-20}}$ ${1.8600812734259762\times 10^{1}}$ ${1.0000000000000001\times 10^{-20}}$ ${1.8600812734259762\times 10^{1}}$
 ${h}$ ${f'(1.5)}$, Central Difference ${1.0000000000000000\times 10^{-2}}$ ${1.8602018344501879\times 10^{1}}$ ${1.0000000000000000\times 10^{-3}}$ ${1.8600824790340198\times 10^{1}}$ ${1.0000000000000000\times 10^{-4}}$ ${1.8600812854816517\times 10^{1}}$ ${9.9999999999999991\times 10^{-6}}$ ${1.8600812735591887\times 10^{1}}$ ${9.9999999999999995\times 10^{-7}}$ ${1.8600812732749716\times 10^{1}}$ ${9.9999999999999995\times 10^{-8}}$ ${1.8600812743407857\times 10^{1}}$ ${1.0000000000000000\times 10^{-8}}$ ${1.8600812623503771\times 10^{1}}$ ${1.0000000000000001\times 10^{-9}}$ ${1.8600814222224926\times 10^{1}}$ ${1.0000000000000000\times 10^{-10}}$ ${1.8600814222224926\times 10^{1}}$ ${1.0000000000000001\times 10^{-11}}$ ${1.8600854190253809\times 10^{1}}$ ${9.9999999999999998\times 10^{-13}}$ ${1.8602897000619123\times 10^{1}}$ ${1.0000000000000000\times 10^{-13}}$ ${1.8589574324323621\times 10^{1}}$ ${1.0000000000000000\times 10^{-14}}$ ${1.8562928971732617\times 10^{1}}$ ${1.0000000000000001\times 10^{-15}}$ ${2.0428103653102880\times 10^{1}}$ ${9.9999999999999998\times 10^{-17}}$ ${0.0000000000000000\times 10^{0}}$ ${9.9999999999999992\times 10^{-18}}$ ${0.0000000000000000\times 10^{0}}$ ${1.0000000000000001\times 10^{-18}}$ ${0.0000000000000000\times 10^{0}}$ ${9.9999999999999998\times 10^{-20}}$ ${0.0000000000000000\times 10^{0}}$ ${1.0000000000000001\times 10^{-20}}$ ${0.0000000000000000\times 10^{0}}$

The first table demonstrates that the complex step approximation improves with decreasing ${h}$. On the other hand, the second table shows that the central difference approximation deteriorates after ${h}$ becomes less than, roughly, ${10^{-12}}$ and finally succumbs to the catastrophic subtractive cancellation.

There is more to be said about using complex variables in numerical differentiation, so stay tuned….