--- title: 'pysersic: Bayesian Inference for Galaxy Structure' url: https://www.emergentmind.com/topics/pysersic type: topic --- # pysersic: Bayesian Inference for Galaxy Structure 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 [2306.05454]. 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 [2412.06957][2508.20266][2603.00380][2509.21236]. ## 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 [2306.05454]. 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 [2510.06315][2509.25335][2407.21257]. 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 [2306.05454]. 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 [2510.06315][2407.21257][2509.25335]. ## 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 [2306.05454]. 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 [2306.05454]. 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` [2510.06315]. 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 [2603.00380]. 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 [2412.06957]. 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 [2412.06957]. ## 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-\(n\) profiles, \(n \ge 3\), are notoriously difficult to render accurately because the central brightness rises steeply as \(r \to 0\), 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 [2306.05454]. 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 \[ \mathcal{F}_r(k) = 2\pi\int_0^\infty I(R) J_0(kR)R \, dR \] varies smoothly with Sérsic index \(n\) and spatial frequency \(k\), while remaining bounded between 0 and 1. The training set was constructed from 10,000 random \((k,n)\) pairs, with \(n\) sampled uniformly from \(0.4\) to \(6.5\) and \(k\) sampled from an exponential distribution with rate \(1/15\), using numerical Hankel transforms computed with the `hankel` package and Ogata quadrature at \(h=10^{-4}\) and \(N=\pi/h\). Symbolic regression was then performed with `pysr`, using L2 loss, maximum complexity 50, and \(2.5\times 10^6\) iterations, while searching over algebraic operators chosen to preserve differentiability and JAX compatibility [2508.20266]. The final emulator was not selected solely by loss minimization. Instead, candidate expressions were evaluated in a JAX setting that timed image rendering, \(\chi^2\) likelihood evaluation, and gradient computation, and the chosen expression was the fastest one with L2 loss below \(2\times 10^{-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 \(R_{\rm eff}\) and \(n\) below \(0.4\%\) and scatter below \(1.7\%\), while real-galaxy tests on 100 HSC-SSP galaxies found median or mean fractional differences below \(1\%\), with scatter \(2.6\%\) for \(R_{\rm eff}\) and \(4.9\%\) for \(n\) [2508.20266]. The same study identifies important caveats. The emulator was designed and validated primarily for \(0.5 \lesssim n \lesssim 6\); training over \(0.56\). 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 [2508.20266]. 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 \(b_n\), deriving a new form via symbolic regression that is stated to be accurate over \(0.110\); square science cutouts are initialized at \(2''\), enlarged iteratively until the nearest segment edge lies at least \(0.5''\) from the cutout boundary, and capped at \(6''\). Neighboring sources within \(2''\) 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 \(10<\mathrm{SNR}<100\) 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 [2603.00380]. The multiband UNCOVER workflow differs by fitting all six wide NIRCam bands jointly on a common \(0.04''\)/pixel grid. Neighboring sources within \(1''\) 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 [2412.06957]. 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 \(0.262''\)/pixel resolution in \(g\), \(r\), and \(z\). Contaminants are identified with `sep` at a \(1.5\sigma\) threshold, cutouts are cropped to seven times the semi-major axis of the central object, the package’s `autoprior` is used for flux and \(R_{\rm eff}\), and about 15% of galaxies require manual reprocessing to correct masking or cropping failures [2509.25335]. ## 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 [2306.05454]. 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 [2508.20266]. The multiband UNCOVER implementation reports a precision gain rather than a raw runtime benchmark. Comparing the width of the F356W \(r_{\rm eff}\) 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 \(<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 \(R_{\rm eff}\) is \(-5\%\) at target \(\mathrm{SNR}=10\) and \(10\%\) at target \(\mathrm{SNR}=5\); by contrast, Sérsic index recovery becomes poor at low S/N, and the authors therefore restrict Sérsic-index analysis to \(\mathrm{SNR}_{\rm opt}>15\) [2412.06957]. 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, \(R_{\rm eff}\), and \(n\) are reported as little or no bias, typically below \(2\%\), at \(\mathrm{SNR}>10\), and structural measurements are summarized as achieving fractional error below \(10\!-\!20\%\) above \(\mathrm{SNR}>10\). 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 \(R_{\rm eff}\), and a factor of 1.5 for \(n\), whereas ellipticity uncertainties are described as almost perfectly calibrated [2603.00380]. 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}\lesssim 25\), and almost all modeled systems with \(\mathrm{SNR}\lesssim 50\) cannot be differentiated from a point source regardless of intrinsic extent. For two-component Sérsic+PSF fits, the extended-component \(R_{\rm eff}\) is not usefully recovered across the explored size range when \(\mathrm{SNR}\lesssim 50\), while the PSF-light fraction \(f_{\rm PSF}\) is argued to be a more robust physically motivated compactness diagnostic than \(R_{\rm eff}\) alone. Even at high SNR, very compact high-\(n\) profiles remain degenerate because the fit can lower \(n\) and compensate with larger \(R_{\rm eff}\) [2509.21236]. 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 [2510.06315][2407.21257]. Crowding and centroid drift motivate explicit neighbor modeling and centroid-based quality cuts [2603.00380][2412.06957]. PSF mismatch, drizzle-induced covariance, and background treatment can dominate the practical error budget in high-resolution imaging [2407.21257][2509.21236]. The Fourier emulator is empirically reliable only over its trained domain in \((k,n)\), and the initial implementation does not address wrap-around aliasing [2508.20266]. ## 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 [2306.05454]. 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 \(\log R_e\) offset of about 0.03 dex in the main text [2603.00380][2407.21257]. The package’s scientific uses are correspondingly varied. In JWST/JADES plus `geko`, pysersic provides the continuum structural layer needed to compare H\(\alpha\) sizes, stellar continuum sizes, and kinematics at \(z\approx4\!-\!6\) [2510.06315]. In the DESI post-starburst study, it underpins the conclusion that massive post-starburst galaxies at \(1