Automatic differentiation

Symbolica supports efficient automatic differentiation using (hyper)dual numbers. The user can specify the complete shape of the dual, i.e., the variables and the depths of the derivatives. Since an optimized dual number has to be created for the specified shape at compile time, this feature is Rust only. During the construction of the dual, loops will be automatically unrolled and unnecessary work is prevented.

An easiest version of duals can be created using:

create_hyperdual_single_derivative!(SingleDual, 3);

which creates a dual in 3 variables and has a single derivative in each. One could represent this dual as: \[ a + b \epsilon_x + c \epsilon_y + d\epsilon_z \] in the variables \(x,y,z\).

SingleDual<T> is generic over a floating-point type T.

A fine-tuned dual that supports higher-order differentiation can be generated with create_hyperdual_from_components!:

create_hyperdual_from_components!(
    Dual,
    [
        [0, 0, 0],
        [1, 0, 0],
        [0, 1, 0],
        [0, 0, 1],
        [1, 1, 0],
        [1, 0, 1],
        [0, 1, 1],
        [1, 1, 1],
        [2, 0, 0]
    ]
);

This dual represents:

\[ a + b \epsilon_x + c \epsilon_y + d\epsilon_z + e \epsilon_x \epsilon_y + f \epsilon_x \epsilon_z + g \epsilon_y \epsilon_z + h \epsilon_x \epsilon_y \epsilon_z + i \epsilon_x^2 \]

We can do computations with dual numbers by defining new dual variables and using them just like regular floating point types:

use symbolica::{
    create_hyperdual_from_components,
    domains::{float::NumericalFloatLike, rational::Rational},
};

fn main() {
    let x = Dual::<Rational>::new_variable(0, (1, 1).into()); // x = 1 + ep_x
    let y = Dual::new_variable(1, (2, 1).into()); // y = 2 + ep_y
    let z = Dual::new_variable(2, (3, 1).into()); // z = 3 + ep_z

    let t3 = x * y * z;

    println!("{}", t3.inv());
}

This yields: \[ \frac{1}{6}\left(1-\epsilon_x-\frac{1}{2} \epsilon_y-\frac{1}{3}\epsilon_z+\frac{1}{2} \epsilon_x \epsilon_y + \frac{1}{3} \epsilon_x \epsilon_z +\frac{1}{6} \epsilon_y \epsilon_z+ \frac{5}{6} \epsilon_x \epsilon_y \epsilon_z+ \epsilon_x^2\right) \]