💻 computer science

Dual Numbers for Arbitrary Order Automatic Differentiation

This paper introduces DNAOAD, a Fortran-based framework that overcomes the memory and scalability limitations of existing nested dual number approaches by utilizing a direct, non-nested representation to efficiently compute arbitrary-order derivatives for high-performance scientific computing.

Original authors: F. Peñuñuri, K. B. Cantún-Avila, R. Peón-Escalante

Published 2026-02-06

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: F. Peñuñuri, K. B. Cantún-Avila, R. Peón-Escalante

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Big Picture: A New Way to Do Math "On the Fly"

Imagine you are a chef trying to perfect a recipe. You want to know exactly how the taste changes if you add a tiny bit more salt, or a tiny bit more sugar. In the world of computer science and engineering, this is called differentiation—figuring out how a result changes when you tweak the input.

For a long time, computers have had two main ways to do this:

Symbolic Math: Like a super-smart student who writes out the entire math formula on a piece of paper. It's exact, but if the recipe gets too complicated, the paper becomes a giant, unreadable mess.
Finite Differences: Like tasting the soup, adding a pinch of salt, tasting it again, and guessing the difference. It's easy to do, but it's not very precise, and if you do it too many times, you ruin the soup (numerical errors).

Automatic Differentiation (AD) is the "Goldilocks" method. It calculates the exact change without needing a giant formula or a messy guess. It works by running the computer program once, but keeping a "shadow" of the calculation that tracks how every single step changes.

The Problem: The "Russian Doll" Trap

The paper focuses on a specific type of AD called Dual Numbers. Think of a dual number as a standard number with a tiny "shadow" attached to it.

Standard Number: 5
Dual Number: 5 + (tiny shadow)

When you do math with these, the shadow automatically updates to tell you the derivative (the rate of change).

The Issue:
Most existing computer programs that use dual numbers to find higher-order derivatives (like the second, third, or hundredth change) use a method called nested dual numbers.

The Analogy: Imagine trying to find the 10th derivative. The computer has to put a dual number inside another dual number, which is inside another, like a set of 10 Russian nesting dolls.
The Result: As you go deeper (higher orders), the computer has to carry around more and more dolls. The memory usage explodes. It's like trying to carry a backpack that gets heavier every time you take a step. Eventually, the computer runs out of space (memory) or gets so slow it crashes.

The Solution: DNAOAD (The Flat Backpack)

The authors of this paper, working at the Universidad Autónoma de Yucatán, created a new tool called DNAOAD.

Instead of using the "Russian Doll" (nested) approach, they built a flat, direct system.

The Analogy: Imagine instead of nesting dolls, you have a single, long row of boxes.
- Box 1 holds the main number.
- Box 2 holds the first change.
- Box 3 holds the second change.
- Box 100 holds the hundredth change.
Why it's better: You don't have to open one box to get to the next. You just reach directly into Box 100. This means the computer doesn't get overwhelmed by "nesting." It uses much less memory and can handle very high orders of derivatives (like the 50th or 100th change) without crashing.

How It Works (The Secret Sauce)

To make this flat system work, the authors used a mathematical tool called the Faà di Bruno formula.

The Metaphor: If you want to know how a complex machine (a function made of smaller functions) changes, you can't just look at the outside. You have to understand how every gear inside turns.
The formula acts like a master instruction manual that tells the computer exactly how to combine the changes of all the little gears to get the final result, without needing to build a giant, nested structure.

What This Paper Actually Does

It introduces a new software package (DNAOAD): It is written in Fortran, a language famous for high-performance scientific computing.
It avoids recursion: It explicitly avoids the "nested" data structures that cause memory problems in other languages (like Julia or C++).
It handles "Arbitrary Order": You can tell the computer, "I want the 50th derivative," and it will do it efficiently, provided your computer has enough RAM.
It includes a library of tools: The package comes with pre-built functions for common math operations (like sine, cosine, exponentials) that automatically know how to handle these "flat" dual numbers.

Real-World Example in the Paper

The paper shows a simple example where they calculate the derivatives of a function involving sine and logarithms.

They set the computer to look for derivatives up to the 5th order.
They feed it a complex number.
The program instantly spits out the original value plus the 1st, 2nd, 3rd, 4th, and 5th changes, all in one go.

They also mention a harder challenge: nesting a function inside itself 1,000 times. Symbolic math would fail here (the formula would be too huge), and standard guessing would be too inaccurate. Their method handles this by keeping the data structure flat and efficient.

Summary

The paper presents DNAOAD, a new way to teach computers how to calculate complex rates of change. Instead of building a tower of nested boxes (which collapses under its own weight), they built a long, flat shelf. This allows scientists and engineers to calculate extremely high-order derivatives efficiently, making it a powerful tool for high-performance scientific computing.

Technical Summary: Dual Numbers for Arbitrary Order Automatic Differentiation

Problem Statement
Automatic Differentiation (AD) is a critical technique for computing derivatives with machine precision, avoiding the truncation errors of finite differences and the expression swell of symbolic differentiation. While forward-mode AD using dual numbers is well-established for first- and second-order derivatives, practical implementations for derivatives of arbitrary order face significant scalability challenges. Existing high-level approaches, such as those in PyTorch (graph-based) or Julia's ForwardDiff (nested dual numbers), suffer from severe memory growth and computational inefficiency as the derivative order increases. Specifically, recursive or nested data structures used to represent higher-order dual numbers lead to rapid memory exhaustion and stack overflows, limiting their applicability in high-performance scientific computing.

Methodology
The authors propose DNAOAD, a Fortran-based framework that computes derivatives of arbitrary order using a direct, non-nested representation of dual numbers. The core methodology involves:

Direct Representation: Instead of recursively nesting dual numbers (e.g., a dual number of order $n$ is not defined as a dual number of order $n-1$ ), the framework defines a dual number of order $n$ as a single structure containing an array of $n+1$ coefficients ( $a_0, \dots, a_n$ ). This avoids recursive data structures entirely.
Faà di Bruno Formula Implementation: To handle the chain rule for composite functions $f(g(x))$ at arbitrary orders, the authors implement the Faà di Bruno formula. This involves computing partial Bell polynomials ( $B_{n,k}$ ) iteratively using dynamic programming (Algorithm 1) rather than recursive calls, ensuring computational efficiency.
Operator Overloading: The framework overloads standard Fortran operators and intrinsic functions (e.g., sin, exp, *, +) to operate on the custom dualzn derived type. For elementary functions without simple closed-form higher derivatives (e.g., arcsin), the system relies on the automatic propagation of derivatives through basic arithmetic operations (addition, multiplication, square root) using the Leibniz rule and Faà di Bruno formula.
Software Architecture: The package is distributed as both a standalone source and an FPM (Fortran Package Manager) module. It includes modules for precision handling, the core dual number type (dualzn_mod), and differential operators (diff_mod) for computing gradients, Hessians, and Jacobians.

Key Contributions

Non-Nested Architecture: The primary contribution is the elimination of recursive data structures in favor of a flat, array-based representation. This fundamentally alters the memory scaling behavior, allowing for the computation of very high-order derivatives without the memory explosion associated with nested approaches.
Arbitrary Order Support: The framework supports derivatives of any order $n$ , limited only by available computational resources and numerical precision, rather than being restricted to low orders (typically 1 or 2) by implementation constraints.
Fortran Implementation: The work provides a practical, high-performance implementation in Fortran, a language widely used in scientific computing, making high-order AD accessible to legacy and modern Fortran codebases.
Comprehensive Operator Set: The package includes overloaded operators for elementary functions, complex arithmetic, and matrix operations, enabling users to write numerical algorithms directly in the algebra of dual numbers.

Results and Performance Claims
The paper demonstrates the framework's capability through example programs, such as computing derivatives of $f(z) = \sin(z)\log(z^2)$ up to the fifth order and differentiating a function nested 1,000 times. The authors claim that by avoiding recursion:

Memory Usage: The approach significantly reduces memory consumption compared to nested dual number implementations.
Scalability: The method overcomes the key scalability limitations of existing methods, enabling the efficient computation of derivatives of "very high order."
Efficiency: The iterative computation of Bell polynomials and the direct access to coefficients via indexing facilitate efficient execution suitable for high-performance applications.

Significance
The authors position this work as the first practical implementation of arbitrary-order dual numbers that avoids recursive or nested representations while remaining scalable. The significance lies in bridging the gap between the theoretical understanding of higher-order AD and its practical application in resource-constrained, high-performance scientific computing. By providing a tool that can handle high-order derivatives without memory exhaustion, DNAOAD enables new possibilities in optimization, sensitivity analysis, and numerical modeling where higher-order derivative information is required but previously inaccessible due to implementation barriers. The work is presented as a foundational step for integrating high-order AD into the Fortran ecosystem.