xad-rs

Exact automatic differentiation for Rust. Forward and reverse mode, first and second order, with name-keyed gradient readback for quant-finance workloads. No finite-difference error, no symbolic manipulation — just the chain rule, applied to your program as it runs.

Unofficial Rust port of the C++ XAD library. Not affiliated with the upstream project.

Documentation

• Theory — derivations and decision-tree guidance: docs/README.md.
• API — generated rustdoc: docs.rs/xad-rs.
• Examples — runnable, cross-checked against analytic answers: examples/.

How it works

Write your numerical code once against the trait Real. Instantiate it at the call site against the concrete mode that matches your problem shape — f64 for no-AD, Jet1 for one directional derivative, Jet1Vec for a full gradient in one forward pass, AReal + Tape for a gradient over many inputs, Jet2/Jet2Vec for second order.

use xad_rs::prelude::*;

fn quadratic<R: Real>(x: &R) -> R {
    x.clone() * x.clone() + R::from(2.0_f64) * x.clone() + R::from(1.0_f64)
}

let v_passive = quadratic(&3.0_f64);                  // primal, no AD
let v_jet1    = quadratic(&Jet1::new(3.0_f64, 1.0));  // forward, 1st order
let v_jet2    = quadratic(&Jet2::variable(3.0_f64));  // forward, 1st + 2nd
// Reverse mode: build inside an active Tape (see below).

Pick a mode

Reverse mode breaks even with forward around n ≈ 4 inputs. For n ≫ 4 (e.g. a 30-input swap pricer), reverse is dramatically faster.

Real is implemented for f64, AReal<f64>, Jet1<f64>, and Jet2<f64>. Jet1Vec / Jet2Vec lack the impl because the trait's From<f64> requires knowing the tangent length; use them directly when you want a full gradient or Hessian in one pass.

By mode — quick recipes

Reverse (gradient via tape)

use xad_rs::{AReal, Tape, math};

let mut tape = Tape::<f64>::new(true);
tape.activate();

let mut x = AReal::new(3.0);
let mut y = AReal::new(4.0);
AReal::register_input(std::slice::from_mut(&mut x), &mut tape);
AReal::register_input(std::slice::from_mut(&mut y), &mut tape);

// f(x, y) = x^2 * y + sin(x)
let mut f = &(&x * &x) * &y + math::ad::sin(&x);
AReal::register_output(std::slice::from_mut(&mut f), &mut tape);
f.set_adjoint(&mut tape, 1.0);
tape.compute_adjoints();

assert!((x.adjoint(&tape) - (2.0 * 3.0 * 4.0 + 3.0_f64.cos())).abs() < 1e-12);
assert_eq!(y.adjoint(&tape), 9.0);

Forward (full gradient in one pass)

use xad_rs::Jet1Vec;

let (x, y) = (Jet1Vec::variable(3.0, 0, 2), Jet1Vec::variable(4.0, 1, 2));
let f = &(&x * &x) * &y;  // x^2 * y

assert_eq!(f.partial(0), 24.0);  // df/dx = 2xy
assert_eq!(f.partial(1),  9.0);  // df/dy = x^2

Second order (gamma along one direction)

use xad_rs::Jet2;

let x = Jet2::variable(2.0_f64);
let y = x * x * x;                            // x^3
assert_eq!(y.first_derivative(),  12.0);      // 3x^2
assert_eq!(y.second_derivative(), 12.0);      // 6x

For the full dense Hessian in one pass, use Jet2Vec via compute_full_hessian (see examples/hessian.rs and the chapter 04 cost model).

Named (greeks by name)

Forward, declare-and-scope:

use xad_rs::{NamedForwardTape, NamedForwardScope};

let mut ft = NamedForwardTape::<f64>::new();
let spot_h   = ft.declare_jet1vec("spot",   100.0);
let strike_h = ft.declare_jet1vec("strike", 105.0);
let scope: NamedForwardScope<f64> = ft.into_scope();

let spot   = scope.jet1vec(spot_h);
let strike = scope.jet1vec(strike_h);
let ratio  = spot / strike;

assert!((ratio.partial("spot") - 1.0 / 105.0).abs() < 1e-14);

Reverse, with a name-keyed gradient as an IndexMap:

use xad_rs::NamedTape;

let mut tape = NamedTape::new();
let x = tape.input("x", 3.0);
let y = tape.input("y", 4.0);
let _registry = tape.freeze();

let f = &(&x * &x) * &y + x.sin();
let grad = tape.gradient(&f);

assert!((grad["x"] - (2.0 * 3.0 * 4.0 + 3.0_f64.cos())).abs() < 1e-12);
assert!((grad["y"] - 9.0).abs() < 1e-12);

Composite Jacobian / Hessian helpers

use xad_rs::{compute_jacobian_rev, compute_hessian};

// f: R^2 -> R^2, f(x, y) = [x*y, x + y]
let _jac = compute_jacobian_rev(&[3.0, 5.0], |v| {
    vec![&v[0] * &v[1], &v[0] + &v[1]]
});

// g: R^2 -> R, g(x, y) = x^2 * y + y^3
let _hess = compute_hessian(&[2.0, 3.0], |v| {
    let x2 = &v[0] * &v[0];
    let y3 = &v[1] * &v[1] * &v[1];
    x2 * &v[1] + y3
});

For the exact Hessian (machine precision, one forward pass instead of n reverse sweeps), reach for compute_full_hessian — it returns a NamedHessian with the value, gradient, and full ndarray::Array2<f64> Hessian.

Examples

Run with cargo run --release --example <name>.

License

AGPL-3.0-or-later, matching the upstream XAD project. See LICENSE.md.

Acknowledgements

• The C++ XAD library — architectural inspiration and source of the financial examples.
• QuantLibAAD — the XAD-instrumented QuantLib build; reference for the AAD-on-quant-finance patterns the financial examples in this crate are modelled after.
• num-traits, indexmap, and ndarray for the underlying primitives.