# How ExoJAX Enables Gradient-Based Optimization for Atmospheric Retrieval

> Discover how ExoJAX enables gradient-based optimization for atmospheric retrieval by making the full forward model differentiable with JAX. Get exact gradients without manual coding.

- Repository: [Hajime Kawahara/exojax](https://github.com/hajimekawahara/exojax)
- Tags: internals
- Published: 2026-03-03

---

**ExoJAX leverages JAX’s automatic differentiation to make the entire atmospheric forward model—from molecular opacities to radiative transfer—fully differentiable, enabling exact gradient computation for retrieval algorithms without manual derivative coding.**

ExoJAX is an open-source Python package for exoplanet atmospheric spectroscopy built entirely on JAX. By implementing opacity calculators and radiative-transfer solvers using pure JAX primitives, the library transforms atmospheric retrieval into a differentiable optimization problem. This architecture allows researchers to apply **gradient-based optimization** and Hamiltonian Monte Carlo to exoplanet spectroscopy with exact analytical gradients flowing through every component of the forward model.

## JAX-First Opacity Calculation for Differentiable Cross-Sections

ExoJAX achieves gradient-based retrieval by ensuring the opacity calculation—the computational bottleneck in atmospheric modeling—is entirely JAX-compatible. The `OpaPremodit` class in [`src/exojax/opacity/premodit/api.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/opacity/premodit/api.py) serves as the primary engine for this capability.

### The diffmode Parameter and Derivative Orders

The `OpaPremodit` constructor accepts a `diffmode` flag that selects the order of JAX-compatible kernels (zeroth, first, or second derivative). This flag is stored at initialization and propagated to low-level kernels built with `jax.grad` and `jax.vmap`, as seen in lines 59-62 of [`src/exojax/opacity/premodit/api.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/opacity/premodit/api.py). When `diffmode` is enabled, the opacity engine preserves the autograd tape through every interpolation and line-shape calculation.

### xsvector and xsmatrix Methods

The core cross-section evaluation methods `xsvector` and `xsmatrix` (lines 87-112) compute cross-section vectors and matrices using only JAX operations: `jnp` array manipulations, `dynamic_slice`, `scan`, and `overlap_and_add`. Because these functions contain no NumPy operations or Python loops, `jax.grad` can differentiate through them automatically. The class selects memory-aware JIT-compatible kernels (`xsvector_close`, `xsmatrix_close`, `xsvector_stitch`, `xsmatrix_stitch`) at initialization (lines 124-138), allowing JAX to compile the entire opacity evaluation once and reuse it for all gradient evaluations.

## Differentiable Radiative-Transfer Solvers

Gradient-based optimization requires gradients to flow not just through opacities, but through the full radiative-transfer (RT) chain from atmospheric parameters to observed spectra.

### ArtTransPure and the JAX Graph

The transmission radiative-transfer class `ArtTransPure` in [`src/exojax/rt/trans.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/rt/trans.py) inherits from `ArtCommon`, which stores all state as JAX arrays. Its `run` method (lines 75-108) builds the geometric chord matrix, computes chord-integrated optical depth, and calls JIT-compiled integration routines (`rtrun_trans_pureabs_simpson` or `rtrun_trans_pureabs_trapezoid`). All operations remain within the JAX ecosystem, ensuring gradients flow backward through the opacity and temperature-pressure profiles. The analogous `ArtEmisPure` class provides the same capability for emission spectroscopy.

Because the RT code never materializes NumPy arrays or breaks the JAX graph, the autograd tape contains the complete computation chain: atmospheric parameters → temperature/pressure profiles → opacity → optical depth τ(λ) → chord integration → observable spectrum.

## Computing Gradients for Atmospheric Retrieval

The library includes integration tests demonstrating full gradient computation for transmission spectra. Here is a complete workflow showing how to define a differentiable forward model and compute its gradient:

```python
import jax
import jax.numpy as jnp
from exojax.opacity import OpaPremodit
from exojax.rt import ArtTransPure
from exojax.database.api import MdbHitran
from exojax.utils.grids import wavenumber_grid

# 1. Build a wavenumber grid (JAX-compatible)

nu_grid, wav, _ = wavenumber_grid(22900.0, 26000.0, Nx=3000,
                                 unit="AA", xsmode="premodit")

# 2. Initialise the RT and opacity objects

art   = ArtTransPure(pressure_top=1e-15, pressure_btm=1e1, nlayer=100)
mdb   = MdbHitran("CO", nu_grid, gpu_transfer=True)
opa   = OpaPremodit(mdb=mdb, nu_grid=nu_grid,
                    auto_trange=[490.0, 510.0],
                    dit_grid_resolution=1.0)

# 3. Define a differentiable forward model

def model(params):
    mmr_CO, mu_fid, T_fid, grav, rad, RV = params
    Tarr = T_fid * jnp.ones_like(art.pressure)          # isothermal T-P

    mmr = art.constant_profile(mmr_CO)
    mu  = art.gravity_profile(Tarr, mu_fid*jnp.ones_like(art.pressure),
                              rad, grav)
    xs  = opa.xsmatrix(Tarr, art.pressure)             # JAX-compatible

    dtau = art.opacity_profile_xs(xs, mmr, opa.mdb.molmass, mu)
    Rp2 = art.run(dtau, Tarr, mu_fid*jnp.ones_like(art.pressure), rad, grav)
    return jnp.sqrt(Rp2)                               # observable radius

# 4. Build a scalar objective (e.g. χ²) and obtain its gradient

def objective(params):
    resid = observed - model(params)
    return jnp.sum(resid**2)

grad_obj = jax.grad(objective)   # ← automatic differentiation

```

The corresponding test file [`tests/integration/unittests_long/transmission/transmission_grad_test.py`](https://github.com/hajimekawahara/exojax/blob/main/tests/integration/unittests_long/transmission/transmission_grad_test.py) (lines 64-71) verifies that these gradients evaluate without NaNs, confirming that derivatives propagate correctly through the opacity and RT layers.

## Integration with jaxopt for Gradient-Based Retrieval

ExoJAX provides high-level examples coupling differentiable forward models to the `jaxopt` optimization library. The test [`tests/endtoend/jaxopt/optimize_spectrum_JAXopt_test.py`](https://github.com/hajimekawahara/exojax/blob/main/tests/endtoend/jaxopt/optimize_spectrum_JAXopt_test.py) (lines 90-103) demonstrates constructing a `jaxopt.GradientDescent` optimizer and running it to minimize residuals of an emission spectrum. Because the forward model is JAX-traced, `jaxopt` receives exact gradients automatically, enabling convergence with significantly fewer forward evaluations than finite-difference approaches.

Here is a simplified retrieval example using `jaxopt`:

```python
import jaxopt
import jax.numpy as jnp

# Assume 'forward' is defined as in the previous example

# and 'obs' contains observed data

def loss(par):
    model_spec = forward(par)
    return jnp.mean((model_spec - obs)**2)

optimizer = jaxopt.GradientDescent(fun=loss,
                                   stepsize=1e-5,
                                   maxiter=500)
init = jnp.array([5e-5, 1000., 800., 6e9])  # [mmr, T, g, R]

params, state = optimizer.run(init)
print("Optimised parameters:", params)

```

This approach works with any JAX-compatible optimizer, including L-BFGS, Adam (via `optax`), or Hamiltonian Monte Carlo (via `numpyro`), because the ExoJAX forward model exposes a standard differentiable interface.

## Summary

- **JAX-first architecture**: Every ExoJAX component—from `OpaPremodit` opacity calculators to `ArtTransPure` radiative-transfer solvers—is implemented using pure JAX primitives.
- **Automatic differentiation**: The `diffmode` parameter in `OpaPremodit` configures derivative orders, while `xsvector` and `xsmatrix` methods maintain differentiability through cross-section calculations.
- **End-to-end gradients**: The `run` method in `ArtTransPure` preserves the JAX computation graph, allowing `jax.grad` to propagate derivatives from observed spectra back to atmospheric parameters (temperature, pressure, mixing ratios).
- **Optimizer compatibility**: Because the full forward model is a single differentiable function **f(params) → spectrum**, users can plug in any JAX-compatible gradient-based optimizer without writing manual derivatives.
- **Verified accuracy**: Integration tests in [`tests/integration/unittests_long/transmission/transmission_grad_test.py`](https://github.com/hajimekawahara/exojax/blob/main/tests/integration/unittests_long/transmission/transmission_grad_test.py) and [`tests/endtoend/jaxopt/optimize_spectrum_JAXopt_test.py`](https://github.com/hajimekawahara/exojax/blob/main/tests/endtoend/jaxopt/optimize_spectrum_JAXopt_test.py) validate that gradients compute correctly and enable successful optimization.

## Frequently Asked Questions

### What makes ExoJAX differentiable compared to traditional atmospheric retrieval codes?

Traditional atmospheric codes often mix NumPy operations, Python loops, and external C/Fortran libraries that break the computation graph. ExoJAX reimplements all numerical routines—opacity evaluation in [`src/exojax/opacity/premodit/api.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/opacity/premodit/api.py) and radiative transfer in [`src/exojax/rt/trans.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/rt/trans.py)—using JAX primitives like `jnp`, `scan`, and `dynamic_slice`. This ensures the autograd tape remains intact from input parameters to output spectra, enabling automatic differentiation via `jax.grad`.

### Which optimization libraries work with ExoJAX for atmospheric retrieval?

Any JAX-compatible optimization library works with ExoJAX. The repository includes tested examples using `jaxopt` (GradientDescent, L-BFGS) in [`tests/endtoend/jaxopt/optimize_spectrum_JAXopt_test.py`](https://github.com/hajimekawahara/exojax/blob/main/tests/endtoend/jaxopt/optimize_spectrum_JAXopt_test.py), but users can also employ `optax` for stochastic gradient descent, `numpyro` for Hamiltonian Monte Carlo, or custom gradient-descent loops. All receive exact analytical gradients automatically because the ExoJAX forward model is a pure JAX function.

### How does the diffmode parameter affect gradient computation in OpaPremodit?

The `diffmode` flag in `OpaPremodit` (lines 59-62 of [`src/exojax/opacity/premodit/api.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/opacity/premodit/api.py)) selects the order of JAX-compatible kernels. When set to first or second derivative modes, the opacity engine uses kernels built with `jax.grad` and `jax.vmap`, allowing automatic differentiation through the line-by-line cross-section calculations. In zeroth mode, the code prioritizes speed over differentiability for forward-only calculations.

### Can ExoJAX compute gradients for both transmission and emission spectroscopy?

Yes. ExoJAX provides `ArtTransPure` for transmission spectra and `ArtEmisPure` for emission spectra, both inheriting from `ArtCommon` in [`src/exojax/rt/common.py`](https://github.com/hajimekawahara/exojax/blob/main/src/exojax/rt/common.py). Both classes store state as JAX arrays and use pure JAX operations in their `run` methods, ensuring gradients flow correctly through the radiative-transfer solvers regardless of the observation geometry.