pysersic: Bayesian Inference for Galaxy Structure

• pysersic is a Python package for Bayesian inference that fits parametric Sérsic models to galaxy images using PSF-aware forward modeling techniques.
• It integrates JAX and NumPyro to deliver gradient-based posterior sampling, optimization, and stochastic variational inference for efficient large-scale galaxy analysis.
• Recent extensions enable multiband morphology inference and accelerated rendering via symbolic-regression emulators, bolstering robust performance across diverse observational regimes.

pysersic is a Python package for Bayesian inference of galaxy structural parameters from imaging data, centered on forward modeling of parametric light profiles—especially Sérsic models—while accounting for the point-spread function (PSF). It was introduced as a pure-Python framework built on jax, with inference implemented in numpyro, so that just-in-time compilation, automatic differentiation, and execution on CPUs, GPUs, or TPUs can be combined with gradient-based posterior sampling, optimization, and stochastic variational inference (Pasha et al., 2023). Subsequent work has extended pysersic from single-band Sérsic fitting to multiband morphology inference, symbolic-regression-based Fourier rendering, and large survey-scale catalog production, while also clarifying its practical limits in low-SNR or PSF-dominated regimes (Miller et al., 2024, Miller et al., 27 Aug 2025, Zhang et al., 27 Feb 2026, Whalen et al., 25 Sep 2025).

1. Definition and scientific scope

At its core, pysersic addresses the standard extragalactic problem of deriving structural measurements from galaxy images by fitting a parametric model directly to the observed pixels. The motivating quantities explicitly identified in the software paper are total flux, effective radius, and Sérsic index, with the broader goal of supporting photometry, morphology studies, structural evolution analyses, recovery of profiles beyond the per-pixel noise limit, and measurements near the image-resolution limit through explicit PSF-aware forward modeling (Pasha et al., 2023).

The package is designed for both single-source and multi-source fitting, with a high-level API and direct access to the underlying probabilistic model. Later applications make clear that, in operational use, pysersic is routinely employed to infer additional geometric parameters such as ellipticity, position angle, and light centroid, and to fit a flat sky term when needed. Representative uses include JWST/JADES continuum size measurements in rest-frame FUV, NUV, and optical bands, Legacy Survey structural measurements in $g$, $r$, and $z$, and HST/WFC3 single-component structural fitting of compact post-starburst galaxies (Danhaive et al., 7 Oct 2025, Asali et al., 29 Sep 2025, Zhang et al., 2024).

The central profile family is the Sérsic law. In the software paper it is written as

$$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$

with $F_{total}$ the total flux, $R_e$ the half-light or effective radius, and $n$ the Sérsic index (Pasha et al., 2023). In later application papers, pysersic is used specifically for one-component 2D Sérsic fitting, rather than bulge-disk decomposition, when the objective is robust global size measurement or homogeneous cross-sample comparison (Danhaive et al., 7 Oct 2025, Zhang et al., 2024, Asali et al., 29 Sep 2025).

2. Computational architecture and inference model

pysersic is written in pure Python and built on jax, which provides just-in-time compilation, automatic differentiation, seamless execution on CPUs, GPUs, and TPUs, and composable transformations including vectorized operations. The rendering module is intentionally separated from the frontend API and inference machinery, so rendering algorithms can be interchanged and tested without changing the probabilistic framework (Pasha et al., 2023).

Inference is implemented in numpyro. The software paper explicitly identifies three inference modes: posterior sampling, especially Hamiltonian Monte Carlo and the No-U-Turn Sampler; optimization for point estimates; and stochastic variational inference, including variational inference with normalizing flows. The same paper emphasizes that gradient availability through jax makes NUTS/HMC attractive for galaxy fitting because the posterior can be moderately to highly dimensional and correlated, while variational methods can reduce the number of likelihood evaluations for large-scale applications (Pasha et al., 2023).

Later studies show that these capabilities are used in distinct ways depending on scale and science case. In JWST/JADES continuum fitting, the continuum morphology is inferred with pysersic in stochastic variational inference flow mode, and the resulting posterior uncertainties are propagated into priors for a second-stage grism analysis with geko (Danhaive et al., 7 Oct 2025). In the Abell 2744 global morphology catalog, the fiducial posterior for every $\mathrm{SNR}>10$ fit is obtained with SVI-MVN, while sources with $10<\mathrm{SNR}<100$ additionally receive MCMC inference with NUTS using 2 chains, 500 warm-up iterations, and 1000 sampling steps (Zhang et al., 27 Feb 2026). In the multiband UNCOVER implementation, posterior sampling is done with NUTS after a NeuTra-style reparameterization, using 4 chains, 750 warm-up steps, and 1500 sampling steps each (Miller et al., 2024).

A major architectural extension is the multiband module added to pysersic in version 0.1.4 for simultaneous fitting of six JWST/NIRCam bands. In that formulation, wavelength dependence is represented with a B-spline,

$S(\lambda) = \sum_i^N c_i B_i(\lambda),$

applied to effective radius, Sérsic index, and ellipticity. The implementation uses a 2nd-order B-spline with 3 knots spanning the full wavelength range, with centroid and position angle shared across bands and flux independent in each band, for a total of 18 free parameters for the primary galaxy (Miller et al., 2024).

3. Rendering, PSF treatment, and mathematical extensions

A persistent technical issue for pysersic is accurate rendering of high-concentration Sérsic profiles. The software paper notes that high-$r$0 profiles, $r$1, are notoriously difficult to render accurately because the central brightness rises steeply as $r$2, making naïve pixelization inaccurate. The initial release therefore provided three rendering algorithms: traditional real-space rendering with central oversampling, Fourier-space rendering via a Gaussian-series representation following Shajib (2019), and a hybrid real/Fourier Gaussian-mixture method based on Lang (2020) to reduce aliasing in purely Fourier-space rendering (Pasha et al., 2023).

A later development replaced one of the main rendering bottlenecks with a symbolic-regression emulator for the radial Fourier transform of the Sérsic profile. The motivating observation was that, for unit total flux and unit effective radius, the radial Fourier transform

$r$3

varies smoothly with Sérsic index $r$4 and spatial frequency $r$5, while remaining bounded between 0 and 1. The training set was constructed from 10,000 random $r$6 pairs, with $r$7 sampled uniformly from $r$8 to $r$9 and $z$0 sampled from an exponential distribution with rate $z$1, using numerical Hankel transforms computed with the hankel package and Ogata quadrature at $z$2 and $z$3. Symbolic regression was then performed with pysr, using L2 loss, maximum complexity 50, and $z$4 iterations, while searching over algebraic operators chosen to preserve differentiability and JAX compatibility (Miller et al., 27 Aug 2025).

The final emulator was not selected solely by loss minimization. Instead, candidate expressions were evaluated in a JAX setting that timed image rendering, $z$5 likelihood evaluation, and gradient computation, and the chosen expression was the fastest one with L2 loss below $z$6, at complexity 39. Implemented as a new rendering scheme in pysersic, the emulator yields galaxy morphology measurements about 2.5 times faster than the package’s default hybrid method with minimal loss in accuracy. Injection-recovery tests on 500 synthetic Sérsic profiles gave mean fractional differences in recovered $z$7 and $z$8 below $z$9 and scatter below $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$0, while real-galaxy tests on 100 HSC-SSP galaxies found median or mean fractional differences below $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$1, with scatter $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$2 for $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$3 and $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$4 for $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$5 (Miller et al., 27 Aug 2025).

The same study identifies important caveats. The emulator was designed and validated primarily for $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$6; training over $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$7 was attempted, but symbolic regression struggled for $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$8. The initial emulator implementation also lacks the aliasing mitigation built into the hybrid Gaussian-mixture renderer, so extended emission near image edges remains a Fourier-rendering concern (Miller et al., 27 Aug 2025).

Another mathematically specific extension appears in the SAGA analysis, which modified pysersic so it could fit Sérsic indices below the earlier lower bound. That study updated the approximation for the Sérsic coefficient $$I(R) \propto F_{total}\exp\left[-\left(\frac{R}{R_e}\right)^{1/n}-1\right],$$9, deriving a new form via symbolic regression that is stated to be accurate over $F_{total}$0, and reports that this approximation is now the default used in pysersic’s PixelRenderer when deriving $F_{total}$1 for a given model (Asali et al., 29 Sep 2025).

4. Operational workflows in major observational programs

pysersic has been deployed in several distinct workflow regimes: band-by-band single-component fitting, multiband joint fitting, catalog-scale survey processing, and structural-measurement support for downstream spectroscopic inference. The implementations differ mainly in cutout preparation, neighbor handling, PSF construction, and the balance between SVI and full MCMC.

Study	Data regime	pysersic role
(Danhaive et al., 7 Oct 2025)	JWST/JADES imaging at $F_{total}$2	Independent one-component Sérsic fits for continuum sizes and morphology priors for geko
(Zhang et al., 2024)	HST/WFC3 F110W imaging of DESI PSBs	Single-Sérsic structural fitting with drizzle-matched empirical PSFs and residual-based checks
(Zhang et al., 27 Feb 2026)	20-band JWST/NIRCam UNCOVER/MegaScience mosaics	Survey-scale band-by-band morphology catalog with SVI-MVN, selective NUTS, and fit-quality flags
(Miller et al., 2024)	Six-band JWST/NIRCam UNCOVER fitting at $F_{total}$3	Joint multiband Bayesian morphology inference with wavelength-linked parameters
(Asali et al., 29 Sep 2025)	Legacy Survey DR9 $F_{total}$4 imaging	Homogeneous structural pipeline for environmental size-mass comparisons

In the JADES-plus-geko application, pysersic is used for one-component 2D Sérsic fitting in selected NIRCam bands corresponding to rest-frame FUV, NUV, and optical. Each band is fit independently. The fits account for PSF convolution using mosaicked WebbPSF-based model PSFs, include a flat sky background term, and mask contaminant sources rather than co-fitting them. The inferred continuum morphology—especially $F_{total}$5, ellipticity, position angle, and centroid—is then used both for size comparisons and for generating priors in the grism-based kinematic model (Danhaive et al., 7 Oct 2025).

In the HST post-starburst analysis, pysersic is embedded in a preprocessing pipeline that uses sep to identify and deblend sources, converts segmentation maps into masks, masks sources that are either 2 magnitudes dimmer than the target galaxy or visually identified as a point source, and dilates masks by 3 pixels using SciPy routines. The PSF treatment is especially elaborate: the authors use the center PSF of the empirical $F_{total}$6 WFC3 grid, inject it into each pre-drizzled exposure, rerun AstroDrizzle with the original drizzle parameters, and extract a drizzle-matched PSF for fitting. The science model is a single Sérsic component plus a flat sky background, and deviations from a perfect Sérsic law are assessed afterward with a residual-based correction following Szomoru et al. (2012) (Zhang et al., 2024).

The Abell 2744 catalog provides the clearest example of pysersic at scale. Sources are selected from the UNCOVER/MegaScience DR3 photometric catalog with use_phot=1 and $F_{total}$7; square science cutouts are initialized at $F_{total}$8, enlarged iteratively until the nearest segment edge lies at least $F_{total}$9 from the cutout boundary, and capped at $R_e$0. Neighboring sources within $R_e$1 are simultaneously fit with their own single-Sérsic components if they satisfy specified magnitude and SNR criteria, while other contaminants are masked with segmentation-based masks dilated by 3 pixels. The posterior strategy combines SVI-MVN for all selected sources and MCMC/NUTS for the $R_e$2 subset, and the outputs are filtered into use_fit=2, use_fit=1, and use_fit=0 classes using masked-pixel, posterior-consistency, boundary-hitting, convergence, residual, centroid-offset, and flux-consistency criteria (Zhang et al., 27 Feb 2026).

The multiband UNCOVER workflow differs by fitting all six wide NIRCam bands jointly on a common $R_e$3/pixel grid. Neighboring sources within $R_e$4 and brighter than one magnitude fainter than the target in the long-wavelength detection image are modeled explicitly; other cataloged sources are masked, and an additional short-wavelength stack mask is built to catch blue sources absent from the long-wavelength detection catalog. The cutouts are background-subtracted using the median value of unmasked pixels, but no sky background parameter is fitted (Miller et al., 2024).

The SAGA implementation shows pysersic in a ground-based survey setting. The authors download co-added Legacy Survey DR9 imaging, inverse-variance maps, and PSFs at native $R_e$5/pixel resolution in $R_e$6, $R_e$7, and $R_e$8. Contaminants are identified with sep at a $R_e$9 threshold, cutouts are cropped to seven times the semi-major axis of the central object, the package’s autoprior is used for flux and $n$0, and about 15% of galaxies require manual reprocessing to correct masking or cropping failures (Asali et al., 29 Sep 2025).

5. Validation, performance, and known failure modes

The software paper’s general performance claim is that fully Bayesian posterior estimation can be obtained in a few minutes or less on modern laptops, with JAX compilation and gradient-based inference presented as the key enablers (Pasha et al., 2023). Later studies supply more detailed benchmarks.

The clearest end-to-end timing result comes from the Fourier-emulator rendering study. Relative to pysersic’s default hybrid renderer, the emulator reduces median runtimes from 64 s to 23.5 s for MCMC, from 36.5 s to 15.3 s for SVI, and from 3.33 s to 1.45 s for MAP optimization, with timings measured on a 2023 MacBook Pro with M2 Pro and including setup, compilation, and save overheads (Miller et al., 27 Aug 2025).

The multiband UNCOVER implementation reports a precision gain rather than a raw runtime benchmark. Comparing the width of the F356W $n$1 posterior from multiband fitting to that from single-band fitting, the median ratio is 0.67, which the paper describes as about 50% tighter constraints on average. The multiband and single-band radius estimates agree at the $n$2 level, so the benefit is increased precision rather than a large shift in inferred size. In degraded-S/N tests, the median fractional difference in recovered $n$3 is $n$4 at target $n$5 and $n$6 at target $n$7; by contrast, Sérsic index recovery becomes poor at low S/N, and the authors therefore restrict Sérsic-index analysis to $n$8 (Miller et al., 2024).

The Abell 2744 catalog adds injection-recovery tests in F200W, F250M, and F444W using ideal single-Sérsic mock galaxies injected into real empty-sky cutouts. For successful fits, median systematic offsets in recovered flux, $n$9, and $\mathrm{SNR}>10$0 are reported as little or no bias, typically below $\mathrm{SNR}>10$1, at $\mathrm{SNR}>10$2, and structural measurements are summarized as achieving fractional error below $\mathrm{SNR}>10$3 above $\mathrm{SNR}>10$4. The same tests also show that pysersic’s nominal posterior uncertainties are not fully calibrated: the uncertainties are too small by about a factor of 2 for flux, a factor of 2 for $\mathrm{SNR}>10$5, and a factor of 1.5 for $\mathrm{SNR}>10$6, whereas ellipticity uncertainties are described as almost perfectly calibrated (Zhang et al., 27 Feb 2026).

The most restrictive failure-mode study concerns JWST little red dots. In simulated F444W cutouts, pysersic and GALFIT both have difficulty recovering input parameters for $\mathrm{SNR}>10$7, and almost all modeled systems with $\mathrm{SNR}>10$8 cannot be differentiated from a point source regardless of intrinsic extent. For two-component Sérsic+PSF fits, the extended-component $\mathrm{SNR}>10$9 is not usefully recovered across the explored size range when $10<\mathrm{SNR}<100$0, while the PSF-light fraction $10<\mathrm{SNR}<100$1 is argued to be a more robust physically motivated compactness diagnostic than $10<\mathrm{SNR}<100$2 alone. Even at high SNR, very compact high-$10<\mathrm{SNR}<100$3 profiles remain degenerate because the fit can lower $10<\mathrm{SNR}<100$4 and compensate with larger $10<\mathrm{SNR}<100$5 (Whalen et al., 25 Sep 2025).

Other limitations recur across applications. Single-Sérsic model misspecification appears in clumpy high-redshift galaxies and tidally disturbed systems, where residuals remain even when the global effective radius is considered robust (Danhaive et al., 7 Oct 2025, Zhang et al., 2024). Crowding and centroid drift motivate explicit neighbor modeling and centroid-based quality cuts (Zhang et al., 27 Feb 2026, Miller et al., 2024). PSF mismatch, drizzle-induced covariance, and background treatment can dominate the practical error budget in high-resolution imaging (Zhang et al., 2024, Whalen et al., 25 Sep 2025). The Fourier emulator is empirically reliable only over its trained domain in $10<\mathrm{SNR}<100$6, and the initial implementation does not address wrap-around aliasing (Miller et al., 27 Aug 2025).

6. Scientific applications and position in the software ecosystem

The software paper situates pysersic alongside galfit, imfit, profit, galight/lenstronomy, PetroFit, and PyAutoGalaxy, but the subsequent literature distinguishes pysersic primarily by its combination of Bayesian posterior inference, JAX acceleration, and flexible rendering backends (Pasha et al., 2023). Later survey papers explicitly contrast it with GALFIT: the Abell 2744 catalog states that GALFIT typically underestimates uncertainties relative to pysersic, while the HST post-starburst study finds no significant systematic bias between the two codes and reports a median $10<\mathrm{SNR}<100$7 offset of about 0.03 dex in the main text (Zhang et al., 27 Feb 2026, Zhang et al., 2024).

The package’s scientific uses are correspondingly varied. In JWST/JADES plus geko, pysersic provides the continuum structural layer needed to compare H$10<\mathrm{SNR}<100$8 sizes, stellar continuum sizes, and kinematics at $10<\mathrm{SNR}<100$9 (Danhaive et al., 7 Oct 2025). In the DESI post-starburst study, it underpins the conclusion that massive post-starburst galaxies at $S(\lambda) = \sum_i^N c_i B_i(\lambda),$0 are compact, have dense cores, and are round in projection (Zhang et al., 2024). In Abell 2744, it serves as the engine for a 29,608-source global morphology catalog spanning observed $S(\lambda) = \sum_i^N c_i B_i(\lambda),$1 (Zhang et al., 27 Feb 2026). In the multiband UNCOVER analysis, it enables joint six-band morphology inference and the reported optical size–mass relation at $S(\lambda) = \sum_i^N c_i B_i(\lambda),$2 (Miller et al., 2024). In the SAGA environmental study, it provides the homogeneous structural measurements required to argue that low-mass satellite galaxies around Milky Way analogs are systematically larger than isolated galaxies at fixed stellar mass (Asali et al., 29 Sep 2025).

Taken together, these studies show pysersic as a modular Bayesian structural-fitting framework whose distinctive contributions are not a new light-profile family but a modern implementation strategy: PSF-aware forward modeling, differentiable rendering, posterior-based uncertainty quantification, and extensibility across single-band, multiband, and survey-scale workflows. The same literature also makes its practical boundaries explicit: performance depends strongly on rendering choice, PSF fidelity, SNR, crowding, and the adequacy of a single-Sérsic description, and in the most compact JWST regimes the correct interpretation may be an upper limit or a PSF-dominance statistic rather than a precise recovered size (Pasha et al., 2023, Miller et al., 27 Aug 2025, Miller et al., 2024, Zhang et al., 27 Feb 2026, Whalen et al., 25 Sep 2025).