Come analyze HEASARC, IRSA, and MAST data in the cloud! The Fornax Initiative is now welcoming all interested beta users.

The Analytic Gradient (VJP) Pipeline

XSPEC's fit-statistic gradient – the derivative of the fit statistic with respect to each thawed parameter – is used in three places:

by fit with method leven (Levenberg–Marquardt), to compute the Gauss–Newton step;
by fit with method migrad (Minuit2), to drive its line search and quasi-Newton Hessian estimate;
by hmc, for every leap-frog sub-step of every NUTS trajectory.

Historically all three used a finite-difference (FD) gradient, which probes the forward model by perturbing each parameter in turn and re-folding the perturbed photon spectrum through the response. This is robust and trivially parallelisable but expensive ( $\sim 2 N$ forward calls per gradient for $N$ thawed parameters) and is dominated by FD step-size noise when the model is tabulated (apec, nei, nsa, kerrbb, ).

This section documents the analytic vector-Jacobian-product (VJP) pipeline that has progressively replaced the FD path. As of the current release, 315 model components have registered analytic gradients – the entire apec/vapec/vvapec and broadened-NEI plasma family, the major absorbers (tbabs, phabs, wabs, varabs, tbvarabs), the line-profile family (gauss, lorentz, kyrline, laor2, the Fe K tabulated lorentzians), the reflection family (pexrav, pexriv, pexmon), relativistic disks (kerrbb, kerrd, slimbh), NS atmospheres (nsa, nsatmos, nsmax, nsmaxg, nsx, carbatm, hatm), AGN SED templates (agnsed, agnslim, qsosed, optxagn), the ionized / warm absorbers (absori, zxipcf, wndabs, plcabs), DEM variants (cemekl, cevmkl, wdem, gadem, cheb6), the Compton family (compTT, compbb, compST, comptb, compPS, eqpair, nteea, nthcomp, thcomp, simpl), and the major multiplicative convolutions (cflux, clumin, cpflux, zashift, partcov, gsmooth, lsmooth). Twenty-two of those are convolutions, which require a chain-rule pass through the upstream model spectrum – see “Convolutions and the chain rule” below.

The analytic and FD paths are combined per parameter: each parameter whose owning component has a registered gradient gets an analytic Jacobian column, and any remaining parameter is finite-differenced and merged into the same gradient (the per-column “hybrid”). A model therefore no longer needs to be fully gradiented to use the analytic path – supplying gradients for the expensive components is enough. Only a few structural cases demote the whole gradient to FD; see “Fallback rules and the override” below.

The two-step gradient

For an additive model component $f_c(E_i; \bm{\theta}_c)$ producing photon flux in energy bin $i$ , the contribution to the fit statistic gradient with respect to its parameter $\theta_{c,k}$ is

$\displaystyle \frac{\partial S}{\partial \theta_{c,k}} \;=\; \sum_{j} \, \frac{... ..._{ji} \, A_i\, \frac{\partial f_c(E_i; \bm{\theta}_c)}{\partial \theta_{c,k}},$

where $S$ is the fit statistic (statistic), $m_j$ is the predicted count rate in detector channel $j$ , $R_{ji}$ is the response matrix, and $A_i$ is the effective area at bin $i$ . The inner sum, $\partial m_j / \partial \theta_{c,k}$ , is what the analytic-gradient pipeline computes. The outer sum is the easy part: for $\chi^2$ it is the bin residual divided by the bin variance; for Cash and the W / PG variants it is a closed form in the count and predicted-rate variables.

The pipeline therefore splits the work into two independently testable pieces:

Per-component photon-flux Jacobian. Each model component supplies a function that returns $\partial f_c / \partial \bm{\theta}_c$ at the same grid of incoming photon energies on which the component computes itself. These functions live in Xspec/src/XSFunctions/ and are registered by name in gradientFunctionMap.cxx; the section below summarises the three file-layout conventions in use.
Response folding and chain-rule combination. Model::analyticFoldedJacobian walks the model tree, collects each component's photon Jacobian, applies the same response and detector area used by the forward pass, and combines absorbers, convolutions and additive parents via the product rule. The result is the $\partial m_j / \partial \theta_{c,k}$ matrix the outer statistic-derivative needs.

The split means a model author writing a new analytic gradient only needs to compute the photon-side derivative – the same thing they would have computed for a hand-checked FD test – and never has to think about the response.

Where the gradient functions live

The pipeline does not care which source file holds the implementation – gradientFunctionMap.cxx registers a name-to-function-pointer association and the linker resolves it – so over the lifetime of the analytic-gradient program three file-layout conventions have accumulated in Xspec/src/XSFunctions/. All three are equivalent at the pipeline level; the choice is one of code-organisation taste.

Colocated with the forward. For closed-form forwards whose gradient is short and obviously parallel to the forward expression, the gradient is defined immediately below the forward in the same source file. Examples include powerLaw.cxx (forward + 3 gradients), tbabs.cxx (5 gradients), cflux.cxx, zashift.cxx, gaussianLine.cxx, logpar.cxx, polconst.cxx, pollin.cxx, polpow.cxx, cutoffPowerLaw.cxx and the large Aped.cxx that holds the whole apec/vapec/vvapec family and their broadened siblings (47 gradients in one file). Colocation makes the two implementations easier to audit – the gradient code and the expression it differentiates are visible side-by-side.
Dedicated per-model <model>Gradient.cxx. For more substantial gradients whose derivation does not fit inline alongside the forward without doubling the source file's size, the gradient gets its own file. Examples include phabsGradient.cxx, bremssGradient.cxx, edgeGradient.cxx, notchGradient.cxx, cyclabsGradient.cxx, smedgeGradient.cxx, wabsGradient.cxx, pegpwrlwGradient.cxx, partcovGradient.cxx, and so on. About two dozen dedicated single-model gradient files exist.
Family <family>FamilyGradient.cxx. Where a group of related models share most of their derivative structure – abundance variants, redshift variants, broadened variants of the same underlying spectrum – the gradients live together in a single family file driven by an X-macro that instantiates each member's ...Gradient(...) function from a common template. Examples include apecCIERSFamilyGradient.cxx (cie/vcie/vvcie and the b- and resonance-scattering variants), demFamilyGradient.cxx, lineProfileFamilyGradient.cxx, comptonReflectionFamilyGradient.cxx, compactObjectFamilyGradient.cxx (NS atmospheres and relativistic disks), ionizedAbsorberFamilyGradient.cxx and roughly a dozen others.

In all three cases the function signature is the same. For a closed-form forward like powerlaw the gradient is

void powerLawGradient(const RealArray& energyArray,
                      const RealArray& params,
                      int spectrumNumber,
                      vector<RealArray>& jacobian,
                      const string& initString)
{
   // jacobian[k][i] holds d f_i / d params[k]
   // (same size grid as the forward fluxArray[i])
   ...
}

regardless of which of the three file-layout patterns hosts it. The same signature is used for multiplicative components – the jacobian array holds $\partial g_c / \partial \bm{\theta}_c$ rather than $\partial f_c / \partial \bm{\theta}_c$ , and Model::analyticFoldedJacobian treats the two cases identically.

For tabulated forwards like apec, the gradient interpolates the same per-element line list the forward uses but with respect to the table's derivative columns; for variable-abundance plasma like vapec the partial derivative with respect to the abundance of element $Z$ is the line-list slice of that element alone (a quantity that drops out of the tabulated forward but is easy to extract on the gradient side). For kerrbb, laor2 and the NS atmospheres – all of which read a 1-D pre-computed file – a five-point central FD against the table's parameter axis returns a derivative that is many orders of magnitude tighter than the FD-on-folded-spectrum path can deliver.

Convolutions and the chain rule

Convolutions are the only category of component that needs the upstream photon spectrum to compute its own gradient. For example cflux multiplies its input by a (parameter-dependent) scalar chosen to set the integrated flux over a band, so

$\displaystyle F_\mathrm{out}(E_i) \;=\; \alpha(\bm{\theta}_\mathrm{cflux};\, F_\mathrm{in}) \;\cdot\; F_\mathrm{in}(E_i)$

and the chain rule gives two contributions to $\partial F_\mathrm{out} / \partial \theta_{u,k}$ for any parameter $\theta_{u,k}$ of an upstream component $u$ :

an own-parameter contribution along the directions of $\bm{\theta}_\mathrm{cflux}$ (the convolution's own parameters);
an upstream-Jacobian contribution that depends on $\partial F_\mathrm{in} / \partial \theta_{u,k}$ .

The convolution-gradient signature XSCCConvGrad takes both the upstream flux $F_\mathrm{in}$ and the upstream Jacobian $\partial F_\mathrm{in} / \partial \theta_u$ as inputs, and returns $F_\mathrm{out}$ , the convolution's own-parameter Jacobian rows, and the upstream-Jacobian-after-convolution rows together. The tree-walk in Model::analyticFoldedJacobian feeds each convolution the upstream Jacobian its forward sibling needed, then glues the convolution's two contribution categories into the combined Jacobian column.

Twenty-two convolutions are on the analytic path: the flux convolutions cflux / clumin / cglumin / cpflux, the shift convolutions zashift / vashift / vmshift / zmshift, the smoothing convolutions gsmooth / lsmooth, the relativistic blurring family kdblur / kdblur2 / rdblur / kerrconv, the reflection convolutions, the Compton-related simpl / thcomp, and the partial-covering partcov.

The registry and the model tree-walk

Analytic gradients register themselves with XSPEC at startup via the same name-keyed map (XSModelFunction::addGradientFunctionPointer, addVJPFunctionPointer) that registers forwards. When a model command names component powerlaw, XSPEC looks up both powerLaw (forward) and powerLawGradient (gradient). A component whose gradient lookup fails is not fatal: its parameters are deferred to the per-column FD path described below. The one exception is an ungradiented convolution, which the tree-walk cannot fold and which therefore demotes the whole gradient to FD.

At gradient-time, Model::analyticFoldedJacobian walks the model tree depth-first, collecting per-component photon Jacobians in topological order so each convolution receives its already- materialised upstream Jacobian. The walk handles:

nested constant groups and Markdown-style parentheses in the model expression (e.g. constant*phabs*(powerlaw + apec));
shared (linked) parameters across components in the same or different sources – each link contributes to its parent column exactly once via the topological order;
gain parameters per response – excluded from the analytic path (the analytic Jacobian then represents only the model-side contribution; see fallback below).

The folded Jacobian is multiplied bin-by-bin by the response matrix, AREASCAL and effective area for each spectrum to produce $\partial m_j / \partial \theta_{c,k}$ , then summed against $\partial S / \partial m_j$ to yield $\partial S / \partial \theta_{c,k}$ , which is what fit and hmc consume.

Fallback rules and the override

The pipeline no longer demotes the whole gradient to FD just because one component is ungradiented. Model::walkGroup instead returns a set of deferred parameters — the shape parameters of an ungradiented additive or multiplicative component (its norm column stays analytic) and the master of any parameter linked across data groups — and the caller (StatManager for leven, MinuitCalcStat for migrad, HMCSampler for hmc) finite-differences just those columns and merges them with the analytic ones. For leven the analytic $\times$ deferred cross terms in the curvature matrix stay correct because beta/alpha are accumulated once, after both column sets are built; for hmc the deferred columns are finite-differenced directly on the statistic.

The whole gradient reverts to FD only when the model tree cannot be folded at all:

an active convolution component lacks a registered gradient – a convolution is a barrier to the tree-walk and, unlike an additive or multiplicative component, cannot be deferred column-by-column (acn/pileup is always in this category);
the model contains a mixing component (mix, amx) or any other structure that analyticFoldedJacobian cannot fold (for example an mdefine expression whose gradient cannot be inferred symbolically);
a parameter is linked to more than one master (a multi-source link such as = p1 + p2), which the per-column FD deferral cannot represent;
a thawed parameter is a response parameter (gain offset or slope) – the analytic path does not yet flow gradients through the response;
the model-string key DISABLE_ANALYTIC_GRAD is set to a true value (set via xset).

The first two cases are structural: for hmc they are a hard refusal (it asks you to use chain) rather than an FD fallback, since a fully gradiented model otherwise takes the reverse-mode VJP fast path and a partly gradiented one the forward hybrid above.

The last item is the user override. At the XSPEC prompt xset DISABLE_ANALYTIC_GRAD yes forces every fit (and HMC trajectory) onto the FD path; this is useful for reproducing historical FD-based fit results bit-for-bit and for isolating an analytic gradient as a suspect when a fit behaves unexpectedly. For a persistent override place the same xset line in ~/.xspec/xspec.rc.

When LM enters the analytic path the user sees a one-time per-session notice at chatter $\geq 10$ summarising the change in path; if a later model load triggers a fallback condition the notice is re-emitted in the other direction.

The audit harness

A standalone audit program at Xspec/src/tools/audit_gradients cross-checks every registered analytic gradient against a five-point central FD on the photon-side flux. For each model the audit:

builds the forward + analytic gradient at the model's default parameter values (with per-model overrides for the few cases whose defaults sit on a numerical singularity or grid edge);
builds a five-point central FD with adaptive step size on the same parameter values;
reports the column- $L_\infty$ relative error against the FD result, with an absolute floor of $10^{-9} \cdot \vert\mathrm{flux}\vert _\mathrm{peak}$ to avoid false-flagging FP-noise columns.

A relative-error threshold of $10^{-3}$ counts as PASS; anything worse is flagged for hand investigation. The full machine-readable audit table is at Xspec/src/tools/audit_gradients/gradient_audit.md, with a per-row classification of failures (TABULATION-NOISE-FLOOR for models whose FD probe falls below the table's own precision floor, INVESTIGATE for true disagreements, etc.).

For convolutions a second-pass audit synthesises a powerlaw upstream model (with closed-form $F_\mathrm{in}$ and $\partial F_\mathrm{in} / \partial \alpha$ ) and runs the same FD metric against both the own-parameter rows and the upstream-Jacobian-after-convolution row. Both cflux and zashift pass at machine precision ( $1.4 \times 10^{-11}$ worst relative error).

Known limitations

A short list of currently-known limitations of the analytic pipeline:

Zero-norm saddle pinning. When an additive component has norm $\approx 0$ , the analytic gradient correctly reports $\partial F / \partial($ other params of that component $) \approx 0$ – because the component contributes essentially nothing to the prediction. LM auto-freezes those parameters, which then stay at their initial values. The previous FD path had the same behaviour (FD noise on a zero-contribution component is also $\approx 0$ ). This is an LM property, not a gradient quality issue. method migrad uses a different step strategy and sometimes escapes the saddle; users hitting this should try migrad first.
Free response parameters force FD. When any free parameter is a response parameter (gain offset or slope), the whole fit falls back to FD. Wiring response parameters into the analytic tree-walk is a planned follow-up.
User-defined and mdefine models. Custom model components loaded via lmod, initpackage or defined in-session by mdefine fall back to FD unless the author registers a gradient function. See Appendix C for the local-model build conventions; the gradient function signature is XSCC or XSCCConv for C++ and XSF77 for Fortran, matching the forward signature with an extra jacobian output array.
Cost trade-off on very wide responses. Each LM iteration on the analytic path applies the response matrix once to all Jacobian columns (the is the forward call itself); FD applies the response $\sim 2 N$ times. When is small (typically $\leq 10$ ) and the response is large (a 60k channel XRISM/Resolve fit), the analytic per-iteration cost can actually be higher than FD even though gradient quality is much better. In that regime the LM-iteration-count drop usually makes up for it; benchmark before declaring a regression.

Per-model anomalies (rows the audit harness has flagged for INVESTIGATE or TABULATION-NOISE-FLOOR but not yet resolved) are tracked in gradient_audit.md; users encountering unexpected fit behaviour on those models should consult the table and, as a workaround, issue xset DISABLE_ANALYTIC_GRAD yes for that fit.

NASA | GSFC | Sciences and Exploration

HEASARC

High Energy Astrophysics Science Archive Research Center

The Analytic Gradient (VJP) Pipeline