The analyticity requirement

1 The complex-step derivative

The entire solver rests on having exact Jacobians: Newton’s method advances by the derivative of the residual, and a derivative corrupted even slightly corrupts every search direction, so an inexact Jacobian is not a small error but a slow, unreliable solve. This document explains the engine that supplies those derivatives exactly and at negligible cost in accuracy and the contract it demands of every residual in return. The method is not a marginal improvement over finite differences; it removes their defining weakness entirely, delivering machine-precision derivatives for a step so small that the truncation error vanishes.

The one price is a discipline on how residuals are written — they must be complex-analytic — which is enforced by the smoothness contract and checked by a roll-call test.

1.1 The finite-difference dilemma

A derivative approximated by a finite difference must balance two errors that pull in opposite directions as the step h shrinks. A central difference has a truncation error that falls as h^2, which argues for a small step; but it also subtracts two nearly equal numbers, R(x+h) - R(x-h), and that subtraction amplifies floating-point round-off as \varepsilon_{\text{mach}}/h, which argues for a large one. The total error is therefore minimized at an intermediate step and cannot be driven below an irreducible floor, given for the central difference as:

\text{error} \;\sim\; C_1 h^2 + \frac{C_2\,\varepsilon_{\text{mach}}}{h}, \qquad \text{minimized near } h \sim \varepsilon_{\text{mach}}^{1/3},

where \varepsilon_{\text{mach}} \approx 2.2\times10^{-16} is the machine epsilon and C_1, C_2 depend on the function’s derivatives. The best attainable relative accuracy is around 10^{-10} to 10^{-8}, worse for an automated solver, and the optimal h depends on the unknown magnitude of the higher derivatives, so no fixed step is safe across the whole state space.

1.2 The complex-step formula

The complex-step method removes the subtraction entirely by perturbing into the imaginary direction. For a residual R that is analytic in a neighbourhood of the real axis, a Taylor expansion about the real point x with an imaginary step \mathrm{i}h_{\text{cs}} is given as:

R(x + \mathrm{i}h_{\text{cs}}) \;=\; R(x) + \mathrm{i}h_{\text{cs}}\,R'(x) - \tfrac{1}{2}h_{\text{cs}}^2\,R''(x) - \tfrac{\mathrm{i}}{6}h_{\text{cs}}^3\,R'''(x) + \cdots,

where h_{\text{cs}} is a real step size and the primes denote derivatives with respect to x. Taking the imaginary part isolates the first derivative, and dividing by the step gives it directly:

R'(x) \;=\; \frac{\Im\,R(x + \mathrm{i}h_{\text{cs}})}{h_{\text{cs}}} + \mathcal{O}(h_{\text{cs}}^2),

where \Im denotes the imaginary part. The decisive feature is that R(x) appears only in the real part and R'(x) only in the imaginary part, so extracting the derivative involves no subtraction of nearly equal quantities and therefore no cancellation error. The step h_{\text{cs}} can consequently be made as small as one likes — 10^{-200} is routine — until the \mathcal{O}(h_{\text{cs}}^2) truncation term is smaller than machine precision and the derivative is exact to the last bit. The contrast with the finite difference is shown in Figure 1: the finite-difference error descends to a floor and then rises again as round-off takes over, while the complex-step error falls to machine precision and stays there for every smaller step.

Code

h_fd = np.logspace(0, -15, 60)
err_fd = np.array([abs((f(x0 + h) - f(x0 - h)) / (2 * h) - truth) / abs(truth) for h in h_fd])
err_fd = np.maximum(err_fd, 1e-18)

h_cs = np.logspace(0, -300, 120)
err_cs = np.array([abs(f(x0 + 1j * h).imag / h - truth) / abs(truth) for h in h_cs])
err_cs = np.maximum(err_cs, 1e-18)

fig = go.Figure()
fig.add_trace(go.Scatter(x=h_fd, y=err_fd, mode="lines", line=dict(color=COLORWAY[1]),
                         name="central finite difference"))
fig.add_trace(go.Scatter(x=h_cs, y=err_cs, mode="lines", line=dict(color=COLORWAY[0], dash="dash"),
                         name="complex step"))
fig.add_hline(y=2.2e-16, line=dict(color="#888888", dash="dot", width=1))
fig.add_annotation(x=1e-250, y=2.2e-16, text="machine ε", showarrow=False, yshift=10,
                   font=dict(size=11, color="#888888"))
fig.update_layout(xaxis_title="step size  h", yaxis_title="relative error")
fig.update_xaxes(type="log", autorange="reversed")
fig.update_yaxes(type="log")
fig

Figure 1: Relative error of a first derivative against step size h, for the central finite difference and the complex step, on the Squire–Trapp function at x = 1.5. The finite difference (solid) descends as h^2 to a floor near 10^{-8} and then climbs again as subtractive cancellation dominates — a V that no fixed step escapes. The complex step (dashed) has no subtraction, so it falls to machine precision and stays there for every smaller step, down to h = 10^{-300}. The exact derivative is taken as the complex step at h = 10^{-200}.

The formula’s power comes with a precise condition: the Taylor expansion into the imaginary direction is valid only where R is complex-analytic, so the residual must be built from operations that remain analytic in a neighbourhood of the real axis. Any construct that branches on the flow state, e.g. an absolute value, a sign, a minimum or maximum, an if on a solution variable, destroys analyticity: it takes its decision on the real part and discards or corrupts the imaginary seed travelling alongside, returning a derivative that is silently wrong. This is exactly why every residual kernel in the framework is assembled from the regularized primitives of the smoothness contract rather than from raw branching, and why a static detector flags any residual whose complex-step and finite-difference derivatives disagree, as this is the signature of a hidden non-analytic operation (test: test_cs_detector_flags_nonanalytic).

1.3 Validation

Because the derivative is meant to be exact, its validation is stringent: the assembled complex-step Jacobian is checked against a finite-difference Jacobian to tolerance, both at the operating point and away from it, and every residual is verified to be continuously differentiable through the zero-flow state on every edge (tests: test_cs_jacobian_matches_finite_difference, test_cs_jacobian_matches_fd_off_operating_point, test_residual_c1_through_zero_flow_all_edges).

The complex-step engine is thus the reason the Jacobian is never itself a source of error. The discipline it demands, analyticity of every residual, is set out in the smoothness contract; the machinery that assembles the complex-stepped kernels into a sparse operator is assembly.