Faster noise modelling could sharpen studies of exoplanet atmospheres

Measuring an exoplanet’s atmosphere involves looking for small changes in how much starlight the planet blocks at different wavelengths. But neighbouring wavelength channels can share disturbances from the detector or the star itself. Treating those channels as statistically independent can therefore misrepresent the uncertainty in the inferred atmosphere.

A preprint by Mark Fortune, Neale P. Gibson and Benjamin Pope presents faster ways to account for these shared patterns. Implemented in the open-source Python package LUAS, the algorithms make selected two-dimensional statistical models practical for datasets containing millions of measurements. The work was posted on arXiv and submitted to Publications of the Astronomical Society of Australia; it is not presented as a peer-reviewed result.

Modelling patterns rather than removing a single correction

The researchers use Gaussian processes: statistical models that describe how measurements are related without requiring an explicit formula for every disturbance. A mathematical rule called a kernel specifies, for example, how strongly noise at nearby times or wavelengths should be correlated.

For transit spectroscopy, a physical model describes the planet crossing its star, while the Gaussian process describes additional variations. Instead of subtracting one supposedly perfect noise correction, the calculation averages over possible disturbance patterns allowed by the model. This carries uncertainty about those patterns into the estimated planetary parameters.

The obstacle is computational. A general Gaussian process records the relationships between every pair of measurements in a covariance matrix. For one million measurements, storing that matrix in double precision would require about eight terabytes. Standard calculations also grow roughly with the cube of the number of measurements.

Existing shortcuts, including the widely used CELERITE approach, can be much faster for suitable one-dimensional models. Extending those advantages across both wavelength and time has been harder.

Breaking a large calculation into smaller ones

The preprint develops seven related optimisation schemes using Kronecker product algebra, which represents a large covariance matrix through smaller matrices associated with individual dimensions.

One approach, called LUAS Lasrach, transforms the data so that a coupled two-dimensional calculation becomes separate one-dimensional Gaussian-process calculations. Those smaller problems can then use established fast methods and run in parallel.

Another approach targets disturbances that share a time-varying shape across many wavelengths, with different amplitudes in each channel. The authors use matrix transformations to concentrate these shared components into a smaller problem. Multiple components can have different time correlations, without having to satisfy an orthogonality restriction imposed by some earlier methods.

For the supported models, these are exact rearrangements of the likelihood calculation—the probability assigned to the observations given model parameters—not approximations that discard measurements. However, “exact” does not guarantee that the selected noise model accurately describes a real instrument or star.

What the benchmarks show

For a dataset with 1,024 wavelength channels and 16,384 time samples—about 16.8 million measurements—a single-component version of the Multi Low Rank model evaluated the likelihood in 29 milliseconds on an M4 MacBook Air. An earlier approach implemented with TINYGP took about 0.6 seconds in the comparison. The new benchmark also included a structured form of measurement covariance.

These timings are not the duration of a complete atmospheric analysis. They exclude the initial compilation step and the physical transit model, and inference may require hundreds of thousands of likelihood evaluations. Calculating parameter gradients adds further expense. Hardware and model choice matter: the most flexible combinations retain costly cubic calculations and do not achieve the fastest scaling.

For numerical checks, the researchers compared their algorithms with a conventional full-matrix calculation on a smaller grid of 16 wavelengths and 128 time samples. They report relative differences in log-likelihood below ten times double-precision machine precision in these tests, and recommend double-precision arithmetic.

Conditions and next steps

Two restrictions are central. Observations must occupy a complete rectangular grid, although spacing along each axis need not be uniform. The noise amplitudes must also be separable: their variation across the grid must be expressible through separate wavelength-dependent and time-dependent factors, rather than an arbitrary uncertainty for every measurement.

The paper also connects these models to principal component analysis, a common way to identify shared patterns. In a simulated two-signal example, different time correlations helped the Gaussian-process model separate components that PCA mixed together. That is a demonstration under controlled conditions, not proof of reliable separation in all observations.

Applications to JWST spectroscopy, ground-based atmospheric studies and pixel-level modelling remain prospects for further work. The immediate contribution is computational: making a wider range of joint noise models affordable, so astronomers can test whether their atmospheric conclusions remain robust when shared disturbances are properly considered.