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

TL;DR

This paper develops a symbolic regression-based emulator for the Fourier transform of Sersic galaxy profiles, enabling faster and differentiable galaxy profile fitting with minimal accuracy loss.

Contribution

The authors create a computationally efficient, differentiable emulator for the Sersic profile's Fourier transform using symbolic regression, improving fitting speed and scalability.

Findings

Galaxy profile fitting is accelerated by 2.5 times using the emulator.

The Fourier transform varies smoothly with Sersic index and frequency.

The emulator maintains accuracy in both simulated and real galaxy data.

Abstract

Galaxy profile fitting is a ubiquitous technique that provides the backbone for photometric and morphological measurements in modern extragalactic surveys. A recent innovation in profile fitting algorithms is to render, or create, the model profile in Fourier space, which aims to provide faster and more accurate results. However, the most common parameterization, the S\'ersic profile, has no closed form Fourier transform, requiring the use of computationally expensive approximations. In this paper our goal is to develop an emulator to mimic the radial Fourier transform of the S\'ersic profile, for use in profile fitting. We first numerically compute the radial Fourier transform and demonstrate that it varies smoothly as a function of the S\'ersic index and $k$, the spatial frequency coordinate. Using this set of numerically calculated transforms as a training set, we use symbolic regression to discover an equation which approximates its behavior. This ensures the emulator will be based on computationally efficient and differentiable building blocks. We implement this novel rendering method in the pysersic profile fitter, and ensure it is accurate by conducting both injection-recovery tests using model galaxy profiles and applying multiple rendering methods to a real sample of galaxies in HSC-SSP imaging. Crucially, the Fourier emulator rendering technique enables measurements of morphological parameters of galaxies 2.5 times faster than standard methods with minimal loss in accuracy. This increased performance while maintaining accuracy is a step that ensures these tools can continue to scale with the ever-increasing flow of incoming data.

Equations

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

$${\mathcal{F}}_r(k)=2\pi {\int }_0^{\infty }I(R)J_0(kR)RdR$$

$$\widetilde{{\mathcal{F}}_r}(k,n)=\frac{1}{1+e^{G(k,n)}}$$

$$H(k,n)$$

$$J(k,n)$$

$$G(k,n)$$

$$\begin{array}{rl}{a}_0& =2.245374\\ a_1& =0.029371526\\ a_2& =2.1431181\\ a_3& =-3.7275262\\ a_4& =0.091609545\\ a_5& =0.32785136\end{array}$$

Full text

1 Introduction

Galaxy profile fitting describes a forward modeling approach to measuring the structural parameters and photometry of galaxies from their images. In many cases the galaxy is assumed to be axisymmetric, with a parametric form of the radial surface brightness profile. The best fit parameters are used to investigate the shape of the surface brightness profile, size, and total flux of the galaxy. By directly modeling the image, one is able to account for the effect of the point spread function, or PSF. Profile modeling has thus become a ubiquitous tool for modern extragalactic astronomy for investigating the morphological properties of galaxies (e.g. Lange et al., 2015; Mowla et al., 2019; Kawinwanichakij et al., 2021; Ward et al., 2024; Allen et al., 2024). For large surveys it is used to obtain accurate photometry of extended sources (Bernardi et al., 2013). Such modeling is also necessary to accurately model sources in crowded fields and across instruments with vastly different depths and spatial resolutions (Dey et al., 2019; Weaver et al., 2022).

In many of these applications galaxy surface-brightness profiles are parameterized using a Sérsic profile (Sérsic, 1963; Sersic, 1968). This empirical description parameterizes the radial surface brightness profile using the functional form shown below:

[TABLE]

Here, $F_{\rm total}$ is the total flux of the galaxy, $n$ is the Sérsic index and describes the shape of the profile, and $R_{\rm eff}$ is the effective — or half-light — radius, which describes the spatial scale of the galaxy. These three parameters, along with the axis ratio, position angle and central position, are used to parameterize the light profile of galaxies. The introduction of the index, $n$, consolidated the fixed index profiles such as the de Vaucouleurs (de Vaucouleurs, 1948) and exponential into a single generalized parametrization. Given this additional freedom the Sérsic profile is able to provide a reasonable approximation to the surface brightness profiles of almost all galaxies. It is also possible to create non-axisymmetric models from this base profile by adding Fourier modes or warps (Peng et al., 2002; Stone et al., 2023). There are numerous software packages and methods that perform Sérsic profile fitting, each with their own aims and strengths. Some of the most popular examples are: galfit (Peng et al., 2002), imfit (Erwin, 2015), profit (Robotham et al., 2017), astrophot (Stone et al., 2023), pysersic (Pasha and Miller, 2023) and forcepho (Baldwin et al., 2024).

While the Sérsic profile’s ability to characterize real galaxies has made it the most popular parametrization in the literature, the use of Sérsic functions in profile fitting is accompanied by several numerical challenges. The first concern is how to accurately render model images, i.e., how to assign model flux values at each pixel given a set of Sérsic parameters. This is a particular problem for galaxies with large indices, $n\gtrsim 2$. Near the center of these galaxies, the slope of the surface brightness profile becomes very steep, leading to drastic changes on scales smaller than that of a pixel. Naive implementations which assign fluxes based on the midpoint of the pixel will fail for these profiles, as they implicitly assume the flux at the center of the pixel provides a reasonable approximation for the average of that pixel — and for rapidly changing profiles this assumption no longer holds. The most common method to render the centers of profiles accurately is to “oversample.” By evaluating the Sérsic model at a series of sub-pixel locations and averaging them to produce the flux estimate for a given pixel one obtains accurate model images for steep profiles. Oversampling comes at a computational cost as many more calculations need to be performed per pixel rendered, though it can be performed recursively by checking for convergence to better focus computational resources where necessary (Robotham et al., 2017; Stone et al., 2023).

A promising alternative method is to render the galaxy in Fourier space; with this method, similar accuracy in the central region can be achieved without the need for oversampling (Rowe et al., 2015; Lang, 2020). This provides a more natural way to render the galaxy “above the atmosphere,” i.e., the intrinsic galaxy image before including the effects of the PSF and pixelization. Furthermore, a Fourier transform is required to convolve a pixel-rendered image with a PSF regardless; rendering directly in Fourier space negates the need for this initial computation, improving performance. Unfortunately, despite the theoretical benefits of a Fourier representation, the Sérsic function has no general closed-form Fourier transform111At fixed indices such as $n=0.5$, representing a Gaussian, and $n=1$ representing an exponential the Fourier transforms are known analytically.

A popular workaround, originally suggested in Hogg and Lang (2013), is to approximate the Sérsic profile using a mixture of Gaussians model. Gaussians are particularly nice functions to work with, especially when using Fourier transforms, as the transform of a Gaussian is also Gaussian. Gaussian decompositions of Sérsic profiles can either be pre-computed and stored in a lookup table (Rowe et al., 2015; Lang et al., 2016) or calculated on the fly, such as using the algorithm presented in Shajib (2019). Additionally, if the PSF is also modeled using Gaussians, then the convolution can be performed analytically. The Mixture of Gaussians representation has been implemented in packages for Sérsic profile fitting such as tractor (Lang et al., 2016), pysersic (Pasha and Miller, 2023) and forcepho (Baldwin et al., 2024). Unfortunately this solution is not a panacea — in order to accurately reproduce the Sérsic profiles, a large number of Gaussians ($\gtrsim 10$) must be employed. Each Gaussian is rendered independently, thus (again) increasing the computational cost.

A wholly different approach was put forth in Spergel (2010). The authors propose a novel functional form of the radial profile, defined in terms of Bessel functions, that loosely mimics the Sérsic profile. Their key insight was to specifically design a functional form such that the Fourier transform is convenient to work with. While some studies have utilized this profile (e.g. Li et al., 2023; Tan et al., 2024), given the decades of inertia of use of the Sérsic profile, there is yet to be widespread adoption.

In this paper, we tackle this longstanding problem using a novel approach. Our first insight is recognizing that although the Fourier transform of the Sérsic profile is intractable to calculate analytically, it is possible to calculate numerically. Thus one can easily build up a library of curves that could be used to train an emulator to mimic its behavior. Our second insight is to employ symbolic regression to emulate the Fourier transforms of the Sérsic profile. Symbolic regression is a supervised learning task which aims to search a wide space of functional forms to find the best fit to a given set of data (Angelis et al., 2023). This is usually achieved using a genetic programming approach which begins with a set of simple expressions and combines them to generate more complex equations. For this study we will use symbolic regression to mimic the Fourier transform of the Sérsic profile to build an emulator based on simple, computationally efficient and differentiable building blocks. This last point is particularly salient given the rise of gradient based optimization and sampling techniques being implemented for profile fitting (Stone et al., 2023; Pasha and Miller, 2023; Chen et al., 2024).

The paper is organized as follows: First we calculate the radially symmetric Fourier transform numerically in Section 2. We then emulate these curves utilizing the pysr package (Cranmer, 2023) to perform symbolic regression in Section 3. In Section 4 we implement a rendering method based on this emulator in the pysersic package to ensure its accuracy for profile fitting. We summarize our approach and discuss caveats and future directions in Section 5.

2 The Radially Symmetric Fourier Transform of the Sérsic Profile

In profile fitting we are interested in the radially symmetric Fourier transform. In almost all applications the galaxy is assumed to be axi-symmetric, with only a radial dependence. To use the Fourier transform we must first switch to using polar coordinates and integrate out the dependence on $\theta$. The radially symmetric Fourier transform can be written as:

[TABLE]

where $k$ is the spatial frequency coordinate, and $J_{0}$ is the Bessel function of the first kind. In two dimensions this is equivalent to another closely related integral transform, the Hankel transform.

We use the python package hankel (Murray and Poulin, 2019), to numerically calculate the radial Fourier transform. By using the numerical integration method for Bessel functions derived in Ogata (2005), this package is able to quickly and accurately calculate the Hankel transform for any given function. For all of the calculations performed in this paper, we use a Sérsic profile that is normalized with a total luminosity and effective radius equal to unity. For all calculations performed in this study using hankel, the quadrature bin width, $h$, which controls the accuracy of the integrations is set to $10^{-4}$. Correspondingly, the number of quadrature bins, $N$, we set to $\pi/h$ as suggested in the documentation. This level of accuracy is likely unnecessary, as the reported uncertainties are often lower than machine precision, $<10^{-16}$, however all of the transforms needed for this study are able to be calculated within a few minutes.

Examples of the radial Fourier transform of four example Sérsic profiles with a range of indices are displayed in Figure 1. We highlight $\mathcal{F}_{r}$ (k) for these four Sérsic profiles (left) using two different views of the same data; a log-log (center) and linear (right) plot. We observe that $\mathcal{F}_{r}$ changes smoothly with both, $k$ and $n$, which is encouraging for the prospects of emulating its behavior. There is also important limiting behavior to consider when approximating $\mathcal{F}_{r}$. Firstly it is bounded by $[0,1]$. \mathcal{F}_{r}$$(k,n) tends towards [math] as $k$ increases and it tends towards 1 as $k$ approaches 0 for all values of $n$. In order to accurately emulate $\mathcal{F}_{r}$of the Sérsic profile it is crucial to match these bounds and limiting behaviors.

3 Finding an Approximation to $\mathcal{F}_{r}$ $(k)$ using Symbolic Regression

3.1 Fitting Procedure

The goal of symbolic regression is to perform a search of the space of analytical expressions for the expression that best matches a given set of input data (Angelis et al., 2023). The search begins with a set of simple functions and combines them to form more complex and expressive functions. The equations are typically evolved using a “genetic” algorithm wherein the function is continually mutated to search for the equations which best match the data. This is a multi-optimization problem where the goal is to find the best fit to the data while minimizing the complexity of the resulting expression. In this study we perform symbolic regression using the python package pysr (Cranmer, 2023). This package implements an evolutionary algorithm which simultaneously evolves multiple populations of equations at varying complexities. At each step the equations are mutated, i.e. operators are added, removed, or swapped or constants are optimized in order to search for the best representation of the data. The result is a pool of best fitting equations at varying levels of complexity.

Our goal is to find an approximation, $\tilde{\mathcal{F}_{r}}$ $(k,n)$, that mimics the behavior radial Fourier transform of the Sérsic profile, $\mathcal{F}_{r}$ $(k,n)$. To facilitate an easier search and ensure the proper limiting behavior observed in Fig. 1, we will impose the following structure on the resulting equation:

[TABLE]

where $G(k,n)$ is an arbitrary function that will be searched for using symbolic regression. This constraint is imposed using the template functionality in pysr222Introduced in version 1.0. This form is a Sigmoid function with the opposite sign in the exponential to match the limits of the limits of $\mathcal{F}_{r}$and enforces the boundedness observed in Fig. 1. Initially in our testing, we did not enforce any structure but found many of the resulting equations eventually converged to a similar form as presented in Eqn. 3. Enforcing this structure a-priori allows for more efficient convergence that results in more accurate final equations. Another option to enforce this boundedness would instead forgo using the template and perform the inverse sigmoid transform before fitting, i.e. performing symbolic regression on $\log\left(\frac{1}{1-\mathcal{F}_{r}}\right)$. However this, in essence, optimizes for fractional uncertainty on $\mathcal{F}_{r}$, and places too much emphasis on the high and low values of $k$. Once the sigmoid function has been applied, minor differences in this regime do not alter the behavior of the resulting equation much. Using the template approach allows us to optimize for absolute uncertainty on $\mathcal{F}_{r}$ rather than fractional uncertainty.

We construct our training set by numerically randomly sampling 10,000 values of $k,n$ and using the hankel python package to calculate $\mathcal{F}_{r}$ $(k,n)$. In this calculation, the flux and $R_{\rm eff}$ of the profile are set to unity. The resulting equation can be trivially rescaled to any value and this normalization allows us to focus the variation in $\tilde{\mathcal{F}_{r}}$ $(k,n)$ with $k$ and $n$. For $n$ we randomly sample values between 0.4 and 6.5, using a uniform distribution. Our target range for $n$ is $0.5-6$, encompassing the majority of real galaxies (Lange et al., 2015; Mowla et al., 2019) and extend the window slightly for the training set to limit the decrease in accuracy near the edges target range. We sample the set of $k$ values from an exponential distribution with a rate of $1/15$. We find that this roughly mimics the distribution of $k$ values when a 200 pixel by 200 pixel image is Fourier transformed. This is around the typical size of cutout used for galaxy image modeling, which is often between 100-200 pixels. (van der Wel et al., 2012; Mowla et al., 2019).

We allow the following basic binary operators: $+,-,*,$, and the following unary operators: $\exp,\log,{\rm square},{\rm sqrt}$, inverse and cube. Setting the inverse as a unary operator, instead of the binary division operator, allows us to better control the constraints to avoid numerical issues. We set up a series of constraints on these operators that attempt to avoid circularly nested functions, i.e. a sqrt directly inside a square, or a log directly inside an exponential. Additionally we do not allow any $+$ or $-$ operators inside the inverse to avoid division by 0. In our testing we found this was necessary to ensure numerically stable equations and did not meaningfully alter the final accuracy of the equations.

We use a L2, or mean square error, loss function and a maximum equation complexity of 50. Complexity is measured as the total number of components in an equation; each constant, operator and variable contributes 1 complexity. We use the batching procedure within pysr with a batch size of 500. In our testing we find this was necessary to explore more equations more efficiently. We use 150 population and 2,000 cycles per iteration while setting weight_optimize parameter to $10^{-3}$ as suggested in the documentation to ensure the constants are optimized. We set the Parsimony and adaptive parsimony scaling to $10^{-8}$ and 500 respectively. These parameters control how much to punish high complexity solutions and we use a lower scaling value to encourage exploration of high complexity, but more accurate, equations. We perform the symbolic search for $2.5\times 10^{6}$ iterations, with the results summarized in Figure 2.

3.2 Selecting an Equation

For many applications, the goal of symbolic regression is to find a simple or interpretable equation. However, our aim is to find an accurate and computationally efficient equation, with little emphasis on readability or simplicity. To understand the properties of the final population of equations in this context, Figure 2 compares the computational efficiency against the accuracy of the final pool of equations. For accuracy we use the L2, or mean squared error, loss for a validation set of 1,000 $k,n$-$\mathcal{F}_{r}$ pairs that are drawn from the same distribution as is used for training. For computational efficiency we measure the execution time to render an 150 x 150 pixel image for a given set of parameters and calculate the $\chi^{2}$ log likelihood, assuming no covariance between pixels, with respect to a reference image. This is implemented in jax, same as pysersic, which allows for auto-differentiation capabilities. In this timing we also include the calculation of gradients of the log likelihood. This setup attempts to mimic the calculations used for sampling and optimization when these emulators are implemented into pysersic. This timing is repeated 500 times for each equation and we report the median. For comparison we repeat the timing procedure for the hybrid renderer, which serves as the default in pysersic and find an average execution time of 1.66 ms.

In general we find there is an inverse correlation between accuracy and computational time. The more complex equations are generally more accurate but require more calculations and therefore have longer execution times. However even these most complex equations still have an execution time over twice as fast as the hybrid renderer. We find there is a plateau in this relationship. Above a complexity of about 30 the execution time continues to increase without any significant gains in the accuracy. In order to maximize both computational efficiency and accuracy we aim to select a point near the beginning of this plateau. Specifically we select the equation with the lowest execution time that has L2 loss $<2\times 10^{-6}$. This yields the equation with a complexity of 39, which is displayed below, in Equation. 6.

[TABLE]

We note there is nothing special about the parts of the equation denoted by $H$ and $J$, we have simply separated them out for clarity in presentation. The constants, $a_{i}$, in Eqn. 6 are given by:

[TABLE]

This equation is an outlier in the plane of execution time and accuracy that makes it particularly suited for our use case. Compared to equations of similar complexity in the final pool there are very few instances of the $k$ evaluation within exponential or square root evaluations. When rendering a profile the number of values of $k$ that need to be calculated scales with the number of pixels and these are computationally expensive operations. Given there is only one value of $n$ to evaluate for each profile the prevalence of this variable has less of an effect on the overall execution time. For the remainder of this paper we will use this equation used to emulate the radial Fourier transform of the Sérsic profile and denote it $\tilde{\mathcal{F}_{r}}$.

Our chosen equation is visualized and compared to the true radial Fourier transform in Figure 3. The emulated equation is shown at four different indices, $n=0.5,1,2.5$ and $4$, alongside the residuals when compared to the true value. For $n=0.5,1$, we include the analytically calculated radial Fourier transform alongside the numerically calculated one. The selected equation provides a good match to all the profiles across the entire $k$ range shown, with residuals constrained to $<10^{-2}$. The scatter accuracy does depend on $n$, with the values of $n<1$ showcasing larger residuals. The emulated equation does diverge from the truth at high $k$ for these profiles as well but this is most evident in the top panel with logarithmic scaling. While both the chosen equation and the numerically calculated $\mathcal{F}_{r}$ trend to zero at high $k$, the behavior is slightly different. The emulator function declines quicker, at least for values of $n<1$.

4 Application to Profile Fitting

4.1 Injection-Recovery Tests

Our final equation is selected based on fitting to the radial Fourier Transform of the Sérsic profile; however, this does not necessarily translate into accuracy in the recovered properties of interest when fitting a galaxy (e.g., $n$, $r_{e}$). To ensure our chosen equation is accurate in real space, we perform a set of injection recovery tests using a “perfect” Sérsic profile.

To begin, we implement this emulator for rendering in pysersic.333This is currently implemented in the branch sersic_fourier_emulator here and will be included in the next release of v0.2 Given the linearity of the Hankel transform, it is straightforward to scale the emulator \tilde{\mathcal{F}_{r}}$$\penalty 10000\ (k,n), which assumes a flux and radii of unity, to any arbitrary values during the fitting process. Our “true” profile will be rendered using a heavily over-sampled pixel space renderer. The entire image is over-sampled by a factor of 10 and the inner 8 pixels are over-sampled by an additional factor of 20.

We render a sample of 500 “true” Sérsic profiles onto a 200 pixel by 200 pixel image with $r_{\rm eff}$ randomly sampled between 1-10 pix, and $n$ between $0.5-6.0$. We additionally sample random values of the axis ratio uniformly between 0 and 0.9, position angle uniformly between 0 and $\pi$ and central position in each coordinate following a Gaussian centered on the center of the image with a width of 0.5 pixels. The flux of each profile is assumed to be 100. We then find the best fit to each profile using the Fourier Emulator rendering scheme to ensure we achieve consistent results. To help stabilize the fitting process we add a small amount of Gaussian noise with an r.m.s of $10^{-4}$ per pixel.

We note that in our initial optimization there was one source which was an outlier from the distribution, with best fit parameters $>20\%$ different from the input values. Upon further investigation we found that the best-fit position angle was exactly 0, the limit of our prior, when the input was $2.66$. This suggested an issue with the minimization. To confirm this we restarted the optimization using a different random seed and found the best-fit parameters were in much better alignment with the input values. We have chosen to only show this second minimization for this source and omit the initial run.

The results of our injection recovery tests are displayed in Figure 4. We find strong agreement between the input parameters and those recovered using the Sérsic Fourier Emulator to render the profiles. For the central position, $x_{c}$ and $y_{c}$, ellipticity and position angle we find agreement to a level better than $10^{-3}$ The average difference between the input and recovered parameters is $<0.4\%$ and the scatter is $<1.7\%$ for $r_{\rm eff}$ and $n$ and $<0.5\%$ for the total flux.

To further investigate the accuracy of the new rendering scheme we show the residuals of the recovery test for $n$ and $r_{\rm eff}$ as a function of the input values. The $r_{\rm eff}$ residuals do not show much structure, with the distribution around the true value remaining relatively consistent as a function of input radius. The residuals of $n$, do show systematic variation as a function of the input index. Most notably the emulator systematically underestimates the index by up to $5\%$ near $n=1$. For values of $n>3$ the index is instead overestimated by $1\%$. This mirrors the results shown in Fig. 3 which indicated the emulator was less accurate for profiles with $n\lesssim 1$.

4.2 Real Galaxies

Along with the injection-recovery tests we also wish to test how well this rendering technique performs in a realistic scenario. For this we compare to the hybrid method implemented in pysersic, which is the default. However both of these rendering methods are imperfect. Our goal with this test is not to assess absolute accuracy, as with the injection-recovery tests, but ensure we are able to provide consistent results with an established technique.

We select 100 galaxies from the Galaxy Mass And Assembly (GAMA) survey, Data Release 4, at $0.4<z<0.5$ in the G02 field (Driver et al., 2022). This field overlaps with the HSC-Subaru Strategic Program imaging program which we will use to perform the morphological fits. Focusing on the $i$ band, science cutouts, weight maps and PSFs are downloaded from the public data release 2 using the unagi python package444https://github.com/dr-guangtou/unagi/tree/master (Aihara et al., 2019). Our aim is to assess the accuracy and performance gain of the rendering method based on the Sérsic Fourier emulator compared to the default hybrid method in pysersic.

For each galaxy we perform the morphological fits with pysersic twice: using the hybrid renderer and using the emulator function $\tilde{\mathcal{F}_{r}}$ $(k,n)$. hybrid is the default rendering method in pysersic and implements the method proposed in Lang (2020) using a series of Gaussians to represent the Sérsic profile. In pysersic the amplitudes of each component are derived using the method described in Shajib (2019). Following Lang (2020), the majority of the Gaussian components are rendered in Fourier space while three with the largest width are rendered in observed space to address the issue of aliasing. For each renderer we use the same input data, prior and setup, using a student-T loss function, with $\nu=3$ and simultaneously fitting a flat sky background.

We perform inference using two methods. The first is stochastic variational inference (SVI) (Hoffman et al., 2013). This method approximates the posterior distribution using a normalizing flow (De Cao et al., 2020). We also use a minimizer to estimate the maximum-a-posteriori (MAP) for each galaxy. It is generally faster but can be less accurate as it will only ever approximate the posterior. The next is MCMC sampling using the No U-Turn sampler implemented in numpyro (Hoffman et al., 2014; Phan et al., 2019) using 2 chains with $1,000$ warm-up and $1,000$ sampling steps. Before sampling we reparameterize the parameters using a multi-variate normal following Hoffman et al. (2019). For these tests we use the default settings in pysersic. We ensure the resulting chains all have effective sample sizes $>300$ for all parameters and $\hat{r}<1.05$ for each renderer (Vehtari et al., 2021). This removes 8 galaxies. Of the remaining galaxies the median effective sample size is 1148 and all values of $\hat{r}$ are $\leq 1.01$. These tests are all performed on a 2023 Macbook Pro with an M2 Pro processor utilizing the CPU.

The results of profile fitting for three example galaxies is shown in Figure 6. In this figure we can directly compare the results when using the Fourier Emulator and hybrid rendering methods. In general we find very little difference; the best fit models for both appear indistinguishable by eye. The recovered parameters are also well matched between the two methods. There are some differences but these are generally smaller than the uncertainties.

A quantitative comparison of the recovered parameters using the different methods is displayed in Figure 7. Similar to Fig. 4 the recovered half-light radius and Sérsic index from the two different rendering methods are directly compared. There are two inset histograms: The upper left the distribution of fractional differences, denoted as $\Delta_{f}$ and the bottom right showcasing the difference with respect to the combined uncertainty, denoted $\Delta_{\sigma}$. The combined uncertainty is calculated by adding that measured from each method in quadrature, assuming they are independent.

We find good agreement between the two methods with the mean fractional difference, calculated using the biweight location (Beers et al., 1990), to be $<1\%$ for both parameters. The scatter in $\Delta_{F}$, calculated using the biweight scale, is $2.6\%$ for $r_{\rm eff.}$ and $4.9\%$ for $n$. We compare the distribution $\Delta_{\sigma}$ to a standard normal which one would expect to find if there were no systematic differences between the two. For $r_{\rm eff}$ we find that $\Delta_{\sigma}$ largely follows a normal distribution with a slight bias towards negative values. For $n$, the distribution of $\Delta_{\sigma}$ shows outliers at $>5\sigma$. These appear to be concentrated at small $n$, near 1, mirroring the results shown in Figure 3. The distribution of $\Delta_{\sigma}$ for $n$ is centered at 0, with a width that appears to be slightly narrower than 1.

We investigate the difference in inference time between the hybrid and Fourier Emulator rendering methods in Figure 8. For both inference techniques, SVI and MCMC, we see a $2.5\times$ improvement in speed when using the Fourier Emulator rendering method. For MCMC sampling specifically this drops the median computational time from 64 sec down to 23.5 seconds. This timing includes all overheads including setup, compilation and well as saving the results. For the SVI method we find the median time decreases from 36.5 seconds to 15.3 seconds, when comparing the hybrid and Fourier Emulator rendering methods respectively. The ratio of median times of 2.4 is similar to that for MCMC sampling. For estimating the MAP, the median time ratio is 2.4, with a median time for 3.33 seconds when using the hybrid renderer and 1.45 seconds for the Fourier Emulator method described in this work.

5 Summary and Discussion

Galaxy profile fitting underpins much of the science performed with extragalactic surveys. The most common parameterization is the Sérsic profile, which can reasonably represent most galaxies but is difficult to work with analytically. Crucially, there is not a general known closed form Fourier transform. In this paper we have explored using symbolic regression to emulate the radial Fourier transform, or Hankel transform, of the Sérsic profile. Our goal is to implement an accurate, computationally efficient, and differentiable rendering scheme for galaxy profile fitting.

To achieve this goal we numerically calculated \mathcal{F}_{r}$$(k,n) for a set of $k$, the spatial frequency coordinate and $n$ the Sérsic index, values using the hankel python package. (Murray and Poulin, 2019). Analyzing the results of our symbolic regression fitting, performed using pysr (Cranmer, 2023), in the context of profile fitting, we selected an equation that is both accurate and computationally efficient. The results of the symbolic regression fitting are shown in Figure 2 and the selected equation in Equation 6.

To test this new rendering method we perform two sets of tests. The first is a set of injection-recovery tests using a heavily over-sampled pixel rendering method as a “perfect” profile. We find the emulator provides consistent results with the average fractional difference in the inferred effective radius and Sérsic index to be $<0.4\%$ for both parameters with scatter $<1.7\%$. However there is systematic variation in the accuracy of the index, with a $-5\%$ offset at $n\approx 1$ and $2\%$ offset at $n>3$. The second is comparing this rendering method to the hybrid method implemented in pysersic applied to the same set of 100 galaxies in the HSC-SSP imaging. This is meant to test how it compares to established techniques in a real-world scenario. Using pysersic, we fit the same set of galaxies using the Sérsic Fourier Emulator rendering method and the default hybrid method based on approximating the Sérsic profile using a mixture of Gaussians. Similar to the injection-recovery tests we found consistent results when comparing the two rendering methods.

The systematic uncertainties imparted by the rendering methods proposed here are below the typical uncertainty of morphological catalogs. In order to assess the typical uncertainty we download two publicly available catalogs from van der Wel et al. (2012); George et al. (2024). The van der Wel et al. (2012) catalog is based on HST data and shows a median fractional uncertainty of 8.9% in the radius and 21.5% in the index. George et al. (2024) uses data from the CLAUDS and HSC-SSP surveys and find a median fractional uncertainties of 7.5% for radius and 11.0% for the index. These are all much greater than the systematic uncertainty implied by our injection-recovery tests in Fig. 4. On a catalog level, the level of accuracy the emulator provides appears to be sufficient, however specific science goals may require more precision.

Comparing the computational runtime, we find that the emulator rendering method performs the MCMC sampling a factor of 2.5 faster than the hybrid method with minimal loss in accuracy. This means on modern laptops MCMC sampling can be performed in less than 30 seconds per galaxy, which results in effective sample sizes $\gtrsim 1,000$. Alternative methods like variational inference can decrease the time required to $\sim 15$ seconds, and optimization to find a best-fitting point in $<1.5$ seconds. These improvements in computational efficiency are necessary to continue to scale scientific pipelines to upcoming surveys with ever increasing numbers of galaxies.

It is important to note that the equation we present in this work is not unique nor likely fully optimized for the task for profile fitting. The loss function we used is based on L2 distance in Fourier space and it is unclear how this metric translates to the accuracy of the rendered image in real space. While we have shown in Fig. 4 that the equation presented in this work is able to accurately render Sérsic profiles, it may be possible to define a different loss function which is directly related to the accuracy of the (final) rendered image. This would focus the optimization on a more practical benchmark. Additionally we have investigated the computational efficiency of the results equations post-hoc, but it may be possible to fold this information in during the fitting process as an additional value to optimize. Continued effort tweaking the symbolic regression procedure or using an entirely different approach will likely result in a more accurate and/or faster emulator.

Additionally, this approach is not without its limitations. First, it is an empirical approximation, and therefore is sensitive to inaccuracies and only valid for the range of parameters it has been trained on. The range of $n$ values we have used, $0.5<n<6$ encompasses the majority of real galaxies (e.g. Lange et al., 2015; Mowla et al., 2019) but it is common to see a higher $n$ limit, often up to $n=8$, used in profile fitting. Initially we attempted to find an emulator over the entire range of $0.5<n<8$ however we found the symbolic regression procedure had difficulty producing accurate results when using$\mathcal{F}_{r}$ $(k,n)$ to perform profile fitting. There were systematic deviations on the scale of $\sim 10\%$ for both high and low values of $n$, $n>6$ and $n<1$ respectively. To ensure an accurate equation which covers the vast majority of realistic galaxy profiles we therefore limited our training set, and thus the applicability of our emulator, to profiles with $n<6$.

Another is aliasing, or “wrapping around” of extended sources near the edge of the image. This is a generic issue of rendering in Fourier space due to the periodic nature of the inverse FFT. Lang (2020) introduces a method to overcome this issue using a series of Gaussians to represent the Sérsic profile and render the largest components in real space. This avoids rendering most of the extended emission, which is the most susceptible to aliasing, in Fourier space therefore minimizing the effects of aliasing. Another option is to render on a larger grid and trim the edges back to the original shape. This simplistic approach simply removes the edges of the cutout where the aliasing is most prominent. However, doing so adds significant computational cost as larger grids need to be used.

In this initial implementation of the emulator rendering method we have not attempted any correction for this type of spatial aliasing. For pysersic, the focus is on modeling single galaxies that are almost always in the center of the cutout, so aliasing is often less of a concern. It may be possible to use the emulator to render in Fourier space up to a given spatial scale and perform the rest of the rendering in real space, similar to the Lang (2020) approach, but this will need to be implemented and tested.

In this work we have introduced a new method of rendering Sérsic profiles by emulating its radial Fourier transform. We use symbolic regression to derive an empirical approximation shown in Eqn. 6 which is implemented in pysersic fitting code. We show that this new method performs measurements 2.5 times faster than the default rendering method with minimal loss in accuracy. This new type of rendering for galaxy profile fitting will help increase computational efficiency to allow pipelines to better scale to the ever increasing amount of data. Importantly, this method retains the differentiability necessary for modern optimization and Bayesian inference algorithms.

Acknowledgments

TBM would like to thank the organizers of the JAXtronomy meeting at the Flatiron Institute where the idea for this paper originated. TBM was supported by a CIERA Postdoctoral Fellowship.

This work used computing resources provided by Northwestern University and the Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA). This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology.

GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is https://www.gama-survey.org/

The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) is led by the astronomical communities of Japan and Taiwan, and Princeton University. The instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. The survey was made possible by funding contributed by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), (Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.