🔭 astrophysics

Using Symbolic Regression to Emulate the Radial Fourier Transform of the Sérsic profile for Fast, Accurate and Differentiable Galaxy Profile Fitting

This paper introduces a symbolic regression-based emulator for the radial Fourier transform of the Sérsic profile that, when implemented in the pysersic fitter, enables galaxy morphological parameter measurements 2.5 times faster than standard methods while maintaining high accuracy and differentiability.

Original authors: Tim B. Miller, Imad Pasha

Published 2026-04-17

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Tim B. Miller, Imad Pasha

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are an astronomer trying to take a perfect portrait of a galaxy. Galaxies aren't just random blobs of light; they have specific shapes, sizes, and brightness patterns that tell us how they were born and how they've evolved. To understand them, scientists use a mathematical "recipe" called the Sérsic profile to describe these shapes.

However, there's a problem: calculating this recipe is like trying to solve a complex puzzle where the pieces keep changing shape. The standard way to do this is slow, like trying to paint a masterpiece by dipping a tiny brush into the paint one single drop at a time. It's accurate, but it takes forever, especially when you have millions of galaxies to photograph.

The Problem: The "Black Box" of Galaxy Math

The most popular way to speed things up is to do the math in a different dimension, called Fourier space. Think of this like translating a book from English to a secret code. In this code, the math becomes much easier to handle.

But here's the catch: The specific "code" for the Sérsic profile (the galaxy shape) doesn't have a simple translation key. It's a "black box." To use the fast Fourier method, scientists have to use clumsy approximations—like trying to guess the shape of a cloud by stacking square blocks on top of each other. It works okay, but it's computationally expensive and requires a lot of "blocks" (calculations) to get it right.

The Solution: The "AI" Translator

The authors of this paper, Tim Miller and Imad Pasha, decided to build a translator for this black box. They didn't try to solve the math manually (because it's impossible). Instead, they used a technique called Symbolic Regression.

Think of Symbolic Regression as a very smart, creative robot that looks at thousands of examples of the "true" galaxy shape and tries to write a simple sentence (an equation) that describes it perfectly.

The Training: They first calculated the "true" shape of the galaxy in the secret code (Fourier space) using a super-accurate, but very slow, computer method. They did this for thousands of different galaxy shapes.
The Discovery: They fed this data to their robot. The robot started guessing equations, testing them, and evolving them (like natural selection) until it found a specific mathematical formula that mimicked the true shape almost perfectly.
The Result: They found a new equation that acts as a shortcut. It's like finding a secret tunnel through a mountain instead of hiking over the peak.

Why This Matters: Speed Without Sacrifice

The authors plugged this new "shortcut" equation into a popular software tool called pysersic and tested it.

The Analogy: Imagine you are baking a cake. The old method was like measuring every grain of sugar individually. The new method is like having a pre-measured, perfectly accurate sugar scoop that works instantly.
The Speed: The new method is 2.5 times faster than the current standard.
The Accuracy: Despite being faster, it is just as accurate. The "cake" (the galaxy model) tastes exactly the same.

Why We Need This

We are entering an era of "Big Data" in astronomy. New telescopes are taking pictures of billions of galaxies. If we use the old, slow method, it would take centuries to analyze all this data. With this new "emulator," we can process data fast enough to keep up with the flood of new images.

The Catch

Like any shortcut, there are a few limits:

Range: The new equation works perfectly for most galaxies (the "normal" shapes), but it gets a little wobbly if you try to use it for extremely weird, extreme shapes.
Edge Cases: Because the math is done in a "wrapped" space (like a video game map where walking off the right edge brings you to the left), very large galaxies near the edge of the image can sometimes get a little distorted. The authors note this but say it's a minor issue for most standard observations.

The Bottom Line

This paper is about finding a clever, mathematically "discovered" shortcut that allows astronomers to analyze the universe's structure much faster, without losing any detail. It's a tool that ensures our telescopes can keep up with the universe's endless stream of new discoveries.

1. The Problem

Galaxy profile fitting is a fundamental technique in extragalactic astronomy for measuring structural parameters (e.g., effective radius, S´ersic index, total flux) from images. The S´ersic profile is the standard parametric model used for this purpose. However, fitting these profiles presents significant numerical challenges:

Lack of Closed-Form Fourier Transform: The S´ersic function has no general closed-form Fourier transform. Since modern fitting often occurs in Fourier space (to handle Point Spread Functions and avoid pixelization artifacts), this necessitates computationally expensive approximations.
Computational Bottlenecks:
- Oversampling: Accurate rendering of steep profiles (high S´ersic index, $n$ ) in real space requires "oversampling" (calculating flux at sub-pixel locations), which is computationally costly.
- Gaussian Mixtures: A common workaround is approximating the S´ersic profile as a mixture of $\gtrsim 10$ Gaussians. While Gaussians have known Fourier transforms, rendering many components independently remains slow.
- Alternative Models: Other functional forms (e.g., Bessel-based) exist but lack widespread adoption due to the inertia of the S´ersic model.
Scalability: As next-generation surveys (e.g., LSST, Euclid) generate massive datasets, current fitting methods struggle to scale while maintaining the differentiability required for modern gradient-based optimization and Bayesian inference (MCMC).

2. Methodology

The authors propose a novel approach: emulating the radial Fourier transform (Hankel transform) of the S´ersic profile using Symbolic Regression.

Numerical Generation of Ground Truth:
- They numerically computed the radial Fourier transform, $F_r(k, n)$ , for a wide range of S´ersic indices ( $n \in [0.4, 6.5]$ ) and spatial frequencies ( $k$ ) using the hankel Python package.
- They established that the transform varies smoothly with $n$ and $k$ and is bounded between 0 and 1.
Symbolic Regression:
- Using the pysr package (based on genetic programming), they searched the space of analytical expressions to find an equation that approximates $F_r(k, n)$ .
- Constraint: To ensure the emulator respects the physical bounds ( $0 \le F_r \le 1$ ) and limiting behaviors ( $F_r \to 1$ as $k \to 0$ ), they enforced a sigmoid structure: $\tilde{F}_r(k, n) = \frac{1}{1 + e^{G(k,n)}}$ .
- Optimization: The regression minimized the L2 loss (Mean Squared Error) on a validation set while balancing equation complexity against execution time.
Implementation:
- The resulting equation was implemented in the pysersic profile fitter.
- The emulator is fully differentiable, allowing it to be used with gradient-based optimizers and MCMC samplers (via JAX).

3. Key Contributions

Derivation of a Closed-Form Emulator: The paper provides an explicit, differentiable analytical equation (Eq. 6 in the paper) that approximates the radial Fourier transform of the S´ersic profile with high fidelity.
Integration into pysersic: The method is integrated into an existing, widely used fitting package, replacing the default "hybrid" Gaussian mixture renderer.
Performance Benchmarking: The authors rigorously tested the method against both synthetic "perfect" profiles and real observational data from the HSC-SSP survey.

4. Results

Accuracy:
- Injection-Recovery Tests: When fitting synthetic galaxies generated with an oversampled "truth" renderer, the emulator recovered parameters with high precision. The average fractional difference was < 0.4% for effective radius ( $r_{eff}$ ) and S´ersic index ( $n$ ), with scatter < 1.7%.
- Systematic Bias: A slight systematic bias was observed where the emulator underestimates $n$ by up to 5% near $n \approx 1$ and overestimates it by ~1% for $n > 3$ . However, these errors are significantly smaller than typical uncertainties in morphological catalogs (which are often >10%).
- Real Data Comparison: Applied to 100 real galaxies in the HSC-SSP survey, the emulator produced results consistent with the default hybrid method. The median fractional difference between the two methods was < 1% for both $r_{eff}$ and $n$ .
Computational Efficiency:
- The emulator rendering method is 2.5 times faster than the standard hybrid Gaussian mixture method.
- MCMC Sampling: Median time per galaxy dropped from 64 seconds to 23.5 seconds.
- Optimization (MAP): Median time dropped from 3.33 seconds to 1.45 seconds.
- Variational Inference (SVI): Median time dropped from 36.5 seconds to 15.3 seconds.

5. Significance and Limitations

Significance:
- Scalability: The 2.5x speedup is critical for processing the petabyte-scale datasets expected from future surveys (e.g., Rubin Observatory/LSST).
- Differentiability: Unlike some approximation methods, this emulator is fully differentiable, enabling efficient gradient-based inference and Bayesian analysis without sacrificing speed.
- Accuracy vs. Speed Trade-off: It achieves a "sweet spot" where accuracy is preserved (errors are below typical catalog uncertainties) while drastically reducing computational cost.
Limitations:
- Parameter Range: The emulator is trained for $0.5 < n < 6$ . While this covers the vast majority of galaxies, it may struggle with extreme bulges ( $n > 6$ or up to $n=8$ ).
- Aliasing: Like all Fourier-space rendering, it is susceptible to aliasing (wrap-around artifacts) for extended sources near image edges. The authors note that while pysersic typically models central galaxies where this is less critical, a hybrid approach (rendering large components in real space) might be needed for edge cases in the future.
- Empirical Nature: The equation is an empirical approximation derived via symbolic regression, not a fundamental physical derivation.

Conclusion:
This work successfully bridges the gap between the need for analytical tractability in Fourier space and the complexity of the S´ersic profile. By using symbolic regression to create a fast, differentiable emulator, the authors provide a tool that significantly accelerates galaxy profile fitting without compromising the accuracy required for modern extragalactic science.