Automatic Differentiation¶

Purpose¶

This specification defines how LuPNT computes derivatives – state-transition matrices, measurement Jacobians, and force-model sensitivities – by automatic differentiation (AD) rather than by hand-coded partials or finite differences. AD is the mechanism behind every \(\Phi = \partial x_f/ \partial x_0\) in Dynamics Models / Numerical Integration, every \(H = \partial z/\partial x\) in GNSS Measurement Model / Measurement Models, and the linearized dynamics/measurement callbacks consumed by the Estimation Filters. The goal is to make the differentiable-type contract, the Jacobian API, and the seeding conventions explicit enough that new models get exact derivatives “for free” and tests can compare against the analytic expressions.

The library and the differentiable scalar are declared in cpp/lupnt/core/definitions.h.

The Differentiable Scalar `Real`¶

LuPNT builds on the autodiff library in forward mode. The fundamental scalar is

\[\texttt{Real} \equiv \texttt{autodiff::real},\]

a first-order forward-mode dual number carrying a value and one directional derivative, \((\,\texttt{val},\ \texttt{grad}\,)\). Every LuPNT vector and matrix type is an Eigen container over Real (Vec3, Vec6, VecX, Mat3, MatX, …), with plain-double counterparts (Vec3d, MatXd, …) for storage and I/O:

// cpp/lupnt/core/definitions.h
namespace ad = autodiff;
using Real = ad::real;                 // forward-mode, first order
using ad::at;  using ad::wrt;  using ad::jacobian;  using ad::derivative;
template <int size> using Vec = Matrix<Real, size, 1>;   // Eigen over Real
template <int r, int c> using Mat = Matrix<Real, r, c>;

Because the models (dynamics ComputeRates, force kernels in environment/forces, measurement Compute, frame rotations, …) are all written in terms of Real arithmetic, a derivative propagates through any of them without additional code. The numeric value of a Real is recovered with .val(); cast<double>() converts a Real Eigen expression to its double values.

Note

Forward mode is first-order only: it yields Jacobians, not Hessians, and the propagated state-transition matrix is therefore the first-order variational \(\Phi = \partial x_f/\partial x_0\) (see Dynamics Models). Because autodiff::real carries a single derivative direction, a full \(m \times n\) Jacobian costs \(n\) forward sweeps (one per input), which jacobian and JacobianParallel (below) organize column by column.

The Jacobian Contract¶

The single AD entry point is autodiff’s jacobian, which evaluates a vector-valued function and its Jacobian in the same call. For \(y = f(x)\) with \(x \in \mathbb{R}^n\), \(y \in \mathbb{R}^m\),

\[J = \frac{\partial y}{\partial x} \in \mathbb{R}^{m\times n}, \qquad \texttt{jacobian}(f,\ \texttt{wrt}(x),\ \texttt{at}(x),\ y,\ J),\]

where wrt(x) names the differentiation variables (all of x, or a single component x(i), or a parameter block p) and at(...) supplies the full evaluation point (which may include additional non-differentiated arguments such as \(t_0, t_f, u\)). The value \(y\) and the Jacobian \(J\) are returned together.

A recurring idiom is to reset the seeds by copying the input through cast<double>() before seeding it, so that the Jacobian is taken cleanly with respect to the requested variables and no stale derivative directions leak in:

VecX x_tmp = x.cast<double>();          // drop any incoming seeds
jacobian(f, wrt(x_tmp), at(x_tmp), y, J);   // J = df/dx at x

State-Transition Matrix (Dynamics)¶

The default Dynamics::Propagate(x0, t0, tf, u, stm) forms the STM by differentiating the (seed-reset) no-STM propagation with respect to the initial state (cpp/lupnt/dynamics/dynamics.cc):

VecX x0_tmp = x0.cast<double>();
jacobian(func, wrt(x0_tmp), at(x0_tmp), xf_tmp, *stm);   // stm = d(xf)/d(x0)

where func re-propagates the state. This is the \(\Phi(t_f,t_0)\) of Dynamics Models; for NBodyDynamics the STM overload requires AD to be enabled (use_ad_). The joint state/parameter overload PropagateWithParams runs two sweeps to build the two sensitivity blocks (used by batch/EKF estimation of force-model or clock parameters):

jacobian(func, wrt(x0_tmp), at(x0_tmp, p0_tmp), xf_tmp, *stm_state);  // d(xf)/d(x0)
jacobian(func, wrt(p0_tmp), at(x0_tmp, p0_tmp), xf_tmp, *stm_param);  // d(xf)/d(params)

Integrator Sensitivity (column-parallel AD)¶

For the integrator sensitivity matrix, JacobianParallel (cpp/lupnt/numerics/math_utils.cc) differentiates the re-propagation function one input at a time, distributing the columns across cores with OpenMP. Each thread seeds a single component wrt(x0(i)) and runs one forward sweep to produce the corresponding column of \(J\):

#pragma omp parallel for
for (int i = 0; i < n; i++) {
  VecX  xi;  MatXd Ji;
  VecX  x0_tmp = x0.cast<double>();
  jacobian(func, wrt(x0_tmp(i)), at(x0_tmp), xi, Ji);   // one column
  J_vec[i] = Ji;
}

Note

JacobianParallel is exact forward-mode automatic differentiation evaluated column-by-column, not a finite-difference approximation: func propagates in Real arithmetic, so the derivative is carried through the integration analytically. The cost is \(n\) forward sweeps (re-runs of the trajectory), which is why the columns are parallelized.

Filter Dynamics and Measurement Jacobians¶

Plain dynamics/measurement models (returning only the value) are lifted into the filter callbacks – which must also return the linearization \(F\) and \(H\) – by wrapping them in jacobian (cpp/lupnt/numerics/filters/filter_utils.cc):

// GetFilterMeasurementFunction: H = dz/dx
FilterMeasurementFunction f = [f_meas](const State& x, MatXd* H, MatXd* R) -> VecX {
  VecX y;  VecX x_tmp = x.cast<double>();
  jacobian(f_meas, wrt(x_tmp), at(x_tmp, R), y, *H);
  return y;
};

// GetFilterDynamicsFunction: F = d(xf)/d(x0)
FilterDynamicsFunction g = [f_dyn](const State& x, Real t0, Real tf, const State* u, MatXd* F) {
  VecX xf;  VecX x_tmp = x.cast<double>();
  jacobian(f_dyn, wrt(x_tmp), at(x_tmp, t0, tf, u), xf, *F);
  return State(xf, x.GetNames(), x.GetUnits());
};

These are exactly the f_dyn_ / f_meas_ callbacks the Estimation Filters consume (\(F\) in the EKF predict, \(H\) in the update). A measurement model that instead differentiates its own raw VecX model in-place – e.g. IslCrosslinkMeasurement (Measurement Models) – calls jacobian directly and deliberately re-seeds a fresh state so that the incoming derivative seeds are preserved for the SchmidtEKF consider-state update.

Differentiability Requirements for Models¶

For AD to yield finite, correct derivatives, model code must stay differentiable through every evaluated branch:

Write model math in Real (not double); use .val() only for logging, comparisons, or when returning to non-differentiated interfaces.
Guard non-smooth boundaries. At the domain edges of asin/acos the derivative is infinite, so the SRP shadow function clamps its arguments to a strictly interior value (ClampUnit, \(\pm(1-10^{-12})\)) so that the SRP contribution to \(\Phi\) stays finite across a grazing eclipse boundary (see the SRP note in Dynamics Models).
Avoid discontinuous control flow keyed on the differentiated variable where a smooth derivative is expected; jacobian returns the derivative of the branch actually taken.

Model Boundaries¶

AD is forward mode, first order (autodiff::real): Jacobians only, no second derivatives; a full \(m\times n\) Jacobian costs \(n\) forward sweeps.
The cast<double>() seed reset is required before an independent jacobian call; omitting it can mix stale derivative directions into the result.
NBodyDynamics::Propagate(..., stm) requires use_ad_ enabled; models that are not written in Real cannot be differentiated this way.
AD differentiates the code as written – non-smooth clamps, table lookups, or branch selections produce the derivative of the realized path, which is why smoothing helpers such as ClampUnit are used at known singular edges.