Automatic_Differentiation_via_Dual

Automatic Differentiation via Dual Numbers

The best etas search relies on accurate and efficient derivative calculations. This is where Fast Optimization, specifically implemented using dual numbers, serves as a powerful feature. While traditional numerical finite difference methods (like forward and central difference) are available, they have limitations that can hinder optimization. This section contrasts these methods with the Automatic Differentiation (AD) approach.

Limitations of Numerical Finite Difference Methods

The forward difference method and central difference method are common numerical techniques for approximating derivatives as follows.

Forward Difference: f'(x) ≈ (f(x + h) - f(x)) / h, with a first-order error term O(h) relative to the step size h.

Central Difference: f'(x) ≈ (f(x + h) - f(x - h)) / 2h, with a second-order error term O(h²) relative to the step size h.

Both methods introduce truncation errors because they truncate the Taylor series expansion. While the central difference method generally has a smaller truncation error (O(h²)) than the forward difference method (O(h)), it requires an extra function evaluation at f(x - h).

Critically, both are susceptible to subtractive or round-off errors. As h approaches zero, the machine’s finite precision can cause inaccurate estimations of (f(x + h) - f(x)) or (f(x + h) - f(x - h)). For example, if the true difference is on the order of 10^-11, but the machine precision is only 10^-9, the calculated difference will be inaccurate.

These problems are especially pronounced with ODE models, common in PK/PD analysis. ODE solver results inherently have uncertainties. Thus, f(x + h), f(x), and f(x - h) inherit these uncertainties, making numerical differences unreliable, especially for small h. Inaccurate derivative estimations can significantly hinder optimization.

The Fast Optimization option (see “Population modeling run options”) employs dual numbers to implement AD, providing a robust solution for accurate and efficient derivative calculations. A dual number extends real numbers by adding a dual component, e, where e² = 0. A dual number is represented as (a, be), or equivalently, a + be, where a is the real part and b is the dual part.

For automatic differentiation, an input x to a function is promoted to a dual number (x, e). When a function y(x) is evaluated with this dual number, it becomes a dual number of the form (y, dy/dx e). The real part is the function's value y(x), and the dual part is the derivative dy/dx multiplied by e.

During the evaluation of y(x, e), all arithmetic operations are overloaded to handle dual numbers according to rules derived from Taylor series expansion and the property e² = 0. Examples include:

(a, be) + (c, de) = (a + c, (b + d)e)

(a, be) * (c, de) = (ac, (ad + bc)e)

exp((y, dy/dx e)) = (exp(y), exp(y) dy/dx e)

These rules ensure that after the function is evaluated with the dual number input, the derivative dy/dx is automatically computed alongside the function value y(x), residing in the dual part.

Using dual numbers for AD offers substantial advantages:

Accuracy: It avoids truncation and subtractive errors, achieving machine-precision accuracy in derivative calculations.

Efficiency: It computes the derivative simultaneously with the function evaluation in a single pass, often faster than finite difference, especially for functions with many inputs.

Robustness: It eliminates the need to choose a step size h, removing the associated complexities and potential for inaccuracies.

The Fast Optimization option is seamlessly integrated with ODE solvers. Variables, numbers, and operators within the solver are promoted to their dual number counterparts. The initial values of parameters being optimized are also promoted to dual numbers. The ODE solver then uses overloaded mathematical operators to handle dual number arithmetic. As a result, the solver outputs a dual number representation of the solution, directly providing both the solution and its derivatives with respect to the input parameters.