nolan (hyperjet)

Const-generic hyperdual numbers for automatic differentiation in Rust

nolan is a forward-mode automatic differentiation library built around const-generic, stack-allocated hyperdual numbers (jets). It provides Jet1<N> for gradients, Jet2<N, H> for full Hessians, and Jet3<N, H, T> for third-order tensors — all without Vec or Box, all generic over the parameter dimension N. The same function body type-checks against f64 (no derivatives), Jet1 (gradients), Jet2 (Hessians), or Jet3 (third- order tensors). Linear algebra, statistics, grids, and angle-wrapping primitives layered on top compose naturally with any jet order, so an orbit-determination solve and its covariance propagation can share the same code path.

The library was developed for astrodynamics — orbit determination, covariance propagation, close-approach sensitivity analysis — but the core types are domain-agnostic and useful anywhere exact, stack-allocated forward-mode derivatives matter: optimization, robotics inverse kinematics, physics simulation, ML gradient checking, sensitivity in scientific computing.

Installation

The crate is published on crates.io as hyperjet. The repo and internal codename stay nolan.

[dependencies]
hyperjet = "1.15.0"

use hyperjet::jets::Jet1;

Internal Empyrean callers can alias the dep back to nolan so existing use nolan::... source keeps working unchanged:

[dependencies]
nolan = { package = "hyperjet", version = "1.15.0" }

use nolan::jets::Jet1;

Types

Jet1<N>

First-order jet: tracks a value and its gradient with respect to N parameters.

use hyperjet::jets::Jet1;

// Create a variable seeded at parameter index 0
let x = Jet1::<3>::variable(2.0, 0);
let y = Jet1::<3>::variable(3.0, 1);

let f = x * x + y;
assert_eq!(f.value, 7.0);     // 2² + 3
assert_eq!(f.grad[0], 4.0);   // ∂f/∂x = 2x = 4
assert_eq!(f.grad[1], 1.0);   // ∂f/∂y = 1

Jet2<N, H>

Second-order jet: tracks value, gradient, and the full Hessian matrix (stored as a lower-triangular array of size H = N(N+1)/2).

use hyperjet::jets::{Jet2, hess_size};

let x = Jet2::<2, { hess_size(2) }>::variable(2.0, 0);
let y = Jet2::<2, { hess_size(2) }>::variable(3.0, 1);

let f = x * x * y;
// f = x²y, value = 12
// ∂f/∂x = 2xy = 12, ∂f/∂y = x² = 4
// ∂²f/∂x² = 2y = 6, ∂²f/∂x∂y = 2x = 4, ∂²f/∂y² = 0

Jet3<N, H, T>

Third-order jet: tracks value, gradient, Hessian, and the full third-order tensor (stored as a lower-triangular array of size T = N(N+1)(N+2)/6).

use hyperjet::jets::{Jet3, hess_size, tens_size};

let x = Jet3::<1, { hess_size(1) }, { tens_size(1) }>::variable(1.0, 0);
let f = x.powi(4);
// f = x⁴, value = 1
// df/dx = 4x³ = 4
// d²f/dx² = 12x² = 12
// d³f/dx³ = 24x = 24
assert_eq!(f.tens[0], 24.0);

Type Aliases

# use hyperjet::jets::{Jet1, Jet2, Jet3};
type Dual = Jet1<1>;                     // Single-variable first derivative
type HyperDual = Jet2<2, 3>;            // Two-variable second derivatives
type HyperHyperDual = Jet3<2, 3, 4>;    // Two-variable third derivatives

Traits

• Differentiable — Base trait: value(), constant(), variable(). Implemented for f64, Jet1, Jet2, Jet3.
• FirstOrder — First derivatives: grad(i). Implemented for Jet1, Jet2, Jet3.
• SecondOrder — Second derivatives: hess(i, j). Implemented for Jet2, Jet3.
• ThirdOrder — Third derivatives: tens(i, j, k). Implemented for Jet3.
• DifferentiableMath — Transcendental functions: sin, cos, tan, asin, acos, atan, atan2, sinh, cosh, tanh, exp, ln, log, sqrt, powf, powi, abs.
• AutoDiff — Marker combining Differentiable, FirstOrder, and DifferentiableMath (the surface every jet order shares — Jet1 implements it, so it does NOT include SecondOrder/ThirdOrder). Use this as a generic bound for order-agnostic code.

Generic Programming

Write functions once that work with f64 (no derivatives), Jet1 (gradients), Jet2 (Hessians), or Jet3 (third-order tensors):

use hyperjet::traits::AutoDiff;

fn gravity<J: AutoDiff>(x: J, y: J, z: J, mu: f64) -> [J; 3] {
    let r = (x * x + y * y + z * z).sqrt();
    let r3 = r.powi(3);
    // `J * f64` is in the `Differentiable` trait bound; `f64 * J` is not
    // generic over J, so we put the scalar on the right.
    [x * (-mu) / r3, y * (-mu) / r3, z * (-mu) / r3]
}

Convenience API: differentiate

For one-shot derivative computations, the differentiate module hides the Jet seeding + extraction boilerplate:

use hyperjet::differentiate::{differentiate1, differentiate2_6, differentiate3_6};

// First derivatives
let (value, grad) = differentiate1([1.5, 2.0], |[x, y]| x * y + x * x);
// value = 5.25, grad = [5.0, 1.5]

// Second derivatives (6-parameter specialization avoids spelling out H)
let (value, grad, hess) =
    differentiate2_6([1.0, 0.5, 0.1, 0.0, 0.0, 0.0], |[x, y, z, _, _, _]| {
        (x * x + y * y + z * z).sqrt()
    });

// Third derivatives
let (value, grad, hess, tens) = differentiate3_6(state, |xs| compute_miss_distance(xs));

Six variants: differentiate1/differentiate2/differentiate3 for scalar \( f: \mathbb{R}^N \to \mathbb{R} \), plus _vec variants for vector-valued \( f: \mathbb{R}^N \to \mathbb{R}^M \) that return the full Jacobian (and higher-order tensors stacked per output). Specialized differentiate2_6, differentiate2_9, differentiate3_6, differentiate3_9 helpers inline the hessian/tensor sizes for the common 6- and 9-parameter state cases.

The wrapper is essentially overhead-free — the seeding + extraction is in the same ballpark as the compute itself (~0.4 ns overhead on a 25 ns Jet1 gravity evaluation, see benchmark_jets.rs).

Runtime-dispatched differentiate_dyn

When the derivative order must be chosen at runtime — e.g., propagate with Jet2 covariance by default, then escalate to Jet3 for skewness diagnostics when a nonlinearity metric trips at a close-approach event — use the dispatched form:

use hyperjet::differentiate::{
    differentiate_dyn_6, AutoDiffFn, Order, Derivatives,
};
use hyperjet::traits::AutoDiff;

// Write the function body once; it works for Jet1, Jet2, or Jet3.
struct MissDistance;
impl AutoDiffFn<6, 1> for MissDistance {
    fn eval<T: AutoDiff>(&self, xs: [T; 6]) -> [T; 1] {
        let [x, y, z, _vx, _vy, _vz] = xs;
        [(x * x + y * y + z * z).sqrt()]
    }
}

let order = if nonlinearity > threshold { Order::Third } else { Order::Second };
let d = differentiate_dyn_6(order, state, &MissDistance);

// Uniform consumption via accessors; hessians/tensors are Option.
let jac = d.jacobian();
if let Some(hess) = d.hessians() { /* use second-order term */ }
if let Some(tens) = d.tensors() { /* use third-order skewness */ }

The hessian and tensor fields are boxed so the enum stays small regardless of dispatched order — First dispatch costs only ~7 ns more than the flat differentiate1_vec.

Linear Algebra

Stack-allocated generic matrix operations for any N:

use hyperjet::linalg::*;

// Solve Ax = b (Gauss-Jordan with scaled partial pivoting, NR §2.5)
let x = mat_solve(&a, &b).unwrap();

// Cholesky decomposition (symmetric positive-definite)
let l = mat_cholesky(&cov).unwrap();  // A = L L^T

// Inverses, eigen, and decompositions
let inv = mat_inv(&a).unwrap();
let d2 = mahalanobis_distance_squared(&x, &mu, &cov).unwrap();
let (v, lambda) = mat_eigenvector_max(&a, 100, 1e-12);
let (eigs, vecs) = mat_symmetric_eigen(&cov).unwrap();   // Jacobi
let kopp = sym_eigenvalues_3(&cov_3x3);                  // closed-form 3×3
let det = mat_det(&a);
let ld = mat_log_det(&a);
let kappa = condition_number(&a);                         // κ₂(A) = σ_max / σ_min
let sigma_max = mat_largest_singular_value(&a, 100, 1e-12);
let trace_cube = mat_trace_cube(&a, &b);                  // Tr((AB)³)

// Rectangular `f64`-only primitives over const-generic shapes
let c = mat_mul::<2, 3, 4>(&a_2x3, &b_3x4);              // 2×3 · 3×4 → 2×4
let at = mat_transpose::<2, 3>(&a_2x3);
let ata = mat_ata::<2, 3>(&a_2x3);                       // Aᵀ A, 3×3
let frob = mat_frobenius::<3, 4>(&a_3x4);

Scaled partial pivoting (NR §2.5) underpins all mat_solve / mat_inv / matN_solve / matN_inv / mat_det paths via the shared hyperjet::linalg::NOLAN_REL_TOL (1e-14) and NOLAN_MIN_SCALE (1e-150) constants. mat3_inv / mat3_solve use a relative determinant guard |det| < REL_TOL · max_entry³; for marginally-conditioned 3×3 inputs prefer mat_inv::<3> (carries the scaled pivot).

Stack-allocated specialised fast paths for 3×3, 6×6, and 9×9 matrices (mat3_*, mat6_*, mat9_*) — generic over T: DifferentiableMath so they work with f64 and Jet types alike.

Covariance regularization (hyperjet::linalg::regularize)

use hyperjet::linalg::regularize::*;

// Project a possibly-indefinite covariance onto the PSD cone with a
// minimum eigenvalue floor (Higham 1988).
let report = nearest_psd::<6>(&cov, 1e-12).unwrap();
println!("clipped {} eigenvalues, max delta = {}",
    report.n_clipped, report.max_clip_magnitude);

// Tikhonov damping: A + αI.
let damped = tikhonov::<6>(&cov, 1e-9);

// Tikhonov with before/after κ₂ — surfaces under-damping at the call site.
let report = tikhonov_with_report::<6>(&cov, 1e-6);
if report.condition_number_after > 1e10 {
    // α was too small; retry with a larger damping factor.
}

The RegularizationReport<N> and TikhonovReport<N> structs are #[must_use] — callers cannot silently drop the diagnostic fields.

Optimization

Gain-ratio Levenberg-Marquardt driver (optimization::lm) — an accept/reject trust-region-style solver with Nielsen damping and Moré diagonal scaling, derived from Madsen-Nielsen-Tingleff (2004), Nielsen (1999), Moré (1978), and the MINPACK lmder structure:

use hyperjet::optimization::lm::{solve_nlls, LMConfig, NLLSEvaluation};

let solution = solve_nlls(
    |x: &[f64; 3]| {
        // Full evaluation: residuals + Jacobian + cost at x
        Ok::<_, MyError>(NLLSEvaluation { residuals, jacobian, cost })
    },
    |x: &[f64; 3]| {
        // Trial cost only (same objective; rejected trials never
        // pay for a Jacobian)
        Ok::<_, MyError>(cost_at(x))
    },
    [0.0; 3],           // initial guess
    &LMConfig::default(),
    None,                // optional Bayesian prior
)?;

println!("x = {:?}, cost = {}", solution.x, solution.cost);

Stateful problems implement ResidualProblem<N> (or, for assembled normal-equation systems such as Schur-reduced fits, SystemProblem<N> via lm::solve_system):

impl CostProblem<6> for MyProblem {
    type Error = MyError;
    fn evaluate_cost(&mut self, x: &[f64; 6]) -> Result<f64, MyError> { ... }

    // Optional: clamp a proposed step. The driver computes the
    // predicted reduction from the clamped step, and a clamped step
    // can never declare convergence by itself.
    fn constrain_step(&mut self, x: &[f64; 6], delta: &mut [f64; 6]) { ... }

    // Two-phase lifecycle: persistent state (caches, diagnostics)
    // commits ONLY on acceptance, so committed state always matches
    // the retained iterate.
    fn on_step_accepted(&mut self, x: &[f64; 6]) { ... }
    fn on_step_rejected(&mut self, x_trial: &[f64; 6]) { ... }
}
impl ResidualProblem<6> for MyProblem {
    fn evaluate(&mut self, x: &[f64; 6]) -> Result<NLLSEvaluation<6>, MyError> { ... }
}
let solution = lm::solve(&mut problem, x0, &config, prior)?;

Features: gain-ratio accept/reject with Nielsen damping (no initial-damping dead zones — convergence is insensitive to tau by construction), Moré running-max diagonal scaling with equilibrated Cholesky solves (mixed-unit parameter vectors spanning ~28 orders of magnitude), step/cost convergence judged against the undamped Gauss-Newton step, opt-in geodesic acceleration (Transtrum-Sethna second-order steps from a problem-supplied directional second derivative — one single-parameter Jet2 evaluation), driver-owned Bayesian priors entering both the normal equations and the compared costs, an explicit per-axis error taxonomy (covariance failure is carried in the solution as a Result, never fabricated zeros), vouched-system termination (CostProblem::assembly_is_exact — a problem serving approximate assemblies, e.g. from a frozen linearization, declares which systems it vouches for, and the driver never ends a solve on one it does not: a criterion firing on an unvouched system forces a fresh assembly and re-tests, so every termination verdict rests on a real evaluation), and bit-deterministic fixed-order arithmetic with a golden-trace regression test.

Statistics

hyperjet::statistics::distributions ships the scalar-input distribution primitives:

use hyperjet::statistics::{ln_gamma, upper_inc_gamma_reg, chi2_sf,
                         normal_pdf, normal_cdf, normal_sf,
                         normal_cdf_difference};

let p = chi2_sf(reduced_chi2 * dof as f64, dof);  // χ² survival
let z = (x - mu) / sigma;
let prob = normal_cdf(z);      // Φ(z), the lower tail
let exceedance = normal_sf(z); // Q(z) = 1 - Φ(z), the upper tail
let inside = normal_cdf_difference(hi, lo); // Φ(hi) - Φ(lo)

normal_cdf, normal_sf and normal_pdf are accurate to a RELATIVE 1e-14 over the whole range where their result is a normal double, with no cutoff: normal_sf(30.0) is 4.9067e-198, not zero. Past about 37.5σ the result is subnormal and only the spacing of the subnormals is left, which is where that guarantee stops. Ask normal_sf for an upper tail rather than forming 1.0 - normal_cdf(x), which has no significant figures left beyond 8σ and none at all beyond 8.3σ. Ask normal_cdf_difference for the probability of a bracket: written out, normal_cdf(37.0) - normal_cdf(30.0) is zero and the function returns 4.9067e-198.

hyperjet::statistics::multivariate ships N-dimensional Gaussian primitives generic over the state dimension:

use hyperjet::statistics::{
    split_gaussian, sigma_points, sample_statistics,
    sigma_points_scaled, weighted_sample_statistics, SigmaPointScaling,
    univariate_split, MAX_SPLIT_COMPONENTS,
};

// Canonical 2N+1 unscaled Julier-Uhlmann sigma points: sample_statistics
// over the returned set round-trips exactly to (mu, cov).
let points = sigma_points::<6>(&mu, &cov).unwrap();
let (mu_back, cov_back) = sample_statistics::<6>(&points).unwrap();

// Merwe scaled unscented transform with a tunable spread (alpha, beta,
// kappa). Propagate `sp.points` through a map, then reconstruct the
// transformed moments with weighted_sample_statistics; for an affine map
// this recovers (mu, cov) exactly.
let sp = sigma_points_scaled::<6>(&mu, &cov, &SigmaPointScaling::merwe()).unwrap();
let (mu_ut, cov_ut) =
    weighted_sample_statistics::<6>(&sp.points, &sp.weights_mean, &sp.weights_cov).unwrap();

// Gaussian mixture decomposition along a chosen direction, from the
// tabulated Vittaldev-Russell splitting library. Weights are unequal and
// fall away from the centre; the mixture reproduces the parent mean and
// covariance exactly, for 1 <= k <= MAX_SPLIT_COMPONENTS (15). A wider k
// is refused by name rather than approximated. The shared sub-covariance
// is deflated along Σe rather than e, so every component stays PSD
// (Σ_k ⪰ s²Σ) even when the split direction is far from an eigenvector.
let components = split_gaussian::<6>(&mu, &cov, &direction, 3).unwrap();

// The underlying univariate entry: weights, means and the shared
// component width s, all in units of the parent standard deviation.
let entry = univariate_split(3).unwrap();

The library replaces an equal-weight split with means at uniformly spaced multiples of the standard deviation, which matched the first two moments and nothing further. The difference is in the tail. Writing the mixture's mass beyond Δσ along the split direction as a fraction of the parent's:

Under the uniform split the tail got worse as k rose; under the library it improves. No k reproduces the deep tail: every component is narrower than the parent, so the mixture's tail always decays faster in the end.

The table is generated, not transcribed. From the repository, where the generator lives (Cargo.toml's include list keeps examples/ and tools/ out of the published crate):

cargo run --release --example generate_split_library > /tmp/split.rs \
    && mv /tmp/split.rs src/statistics/split_library.rs \
    && cargo fmt

The temporary matters: redirecting straight onto the module truncates it before cargo runs, and the crate then will not compile. So does the cargo fmt, which is what makes the emitted table match the committed file byte for byte. tests/split_library.rs re-solves the optimisation again and checks what is committed, and separately reproduces the seven-component library published in the source paper.

Grids

NumPy-semantics endpoint-inclusive grid generators and a clamping linear interpolator:

use hyperjet::grids::{linspace, logspace, linear_clamped};

let lin = linspace(0.0, 1.0, 11);          // [0.0, 0.1, ..., 1.0]
let log = logspace(1e-3, 1e3, 7);          // log-spaced over 6 decades
let y = linear_clamped(&xs, &ys, query_x); // clamped at boundaries

linear_clamped<T> is generic over T: Copy + Add<Output = T> + Mul<f64, Output = T> — so it works for scalar values and user-defined types that impl those traits.

Angles

use hyperjet::angles::{wrap_pi, wrap_2pi, wrap_180, wrap_360};

let residual = wrap_pi(observed_ra - predicted_ra);  // half-open (-π, π]
let lon = wrap_2pi(atan2_result);                    // [0, 2π)
let dec_deg = wrap_180(d - d_ref);                   // (-180°, 180°]

Version

println!("{}", hyperjet::version());
// Tagged release: "1.0.0"
// Development:    "1.0.1-dev+a3f7b2c"
// Dirty:          "1.0.1-dev+a3f7b2c-dirty"