Galmoss: A package for GPU-accelerated Galaxy Profile Fitting

Abstract

We introduce galmoss, a python-based, torch-powered tool for two-dimensional fitting of galaxy profiles. By seamlessly enabling GPU parallelization, galmoss meets the high computational demands of large-scale galaxy surveys, placing galaxy profile fitting in the LSST-era. It incorporates widely used profiles such as the Sérsic, Exponential disk, Ferrer, King, Gaussian, and Moffat profiles, and allows for the easy integration of more complex models. Tested on 8,289 galaxies from the Sloan Digital Sky Survey (SDSS) g-band with a single NVIDIA A100 GPU, galmoss completed classical Sérsic profile fitting in about 10 minutes. Benchmark tests show that galmoss achieves computational speeds that are 6 $\times$ faster than those of default implementations.

1 Introduction

Galaxies, the cosmic building blocks, comprise diverse stellar components such as bulges, discs, bars, spiral arms, and nuclear star clusters. One of the main drivers of extragalactic astronomy is the study of the structural and morphological properties of galaxies from their photometric images, which has been shown to correlate with the galaxies’ formation and evolutionary paths (van der Wel, 2008; Conselice, 2014; Dimauro et al., 2022).

Several approaches have been developed for galaxy structural and morphological analysis. Classical eyeball morphology classifications, such as the Hubble sequence (Hubble, 1926), are generally descriptive and rely on visual inspection. Non-parametric morphological analysis (Ferrari et al., 2015; Rodriguez-Gomez et al., 2019) employs quantitative metrics to represent attributes such as concentration, asymmetry, and smoothness (Conselice, 2003). These metrics are considered robust because they do not depend on any underlying model. In contrast, the analysis of surface brightness profiles of galaxies involves fitting the light distribution with parametric models to interpret and quantify morphology using specific model parameters (Nantais et al., 2013; Zhuang and Ho, 2022). This profile fitting method has evolved from one-dimensional (de Vaucouleurs, 1958; Sersic, 1968) to two-dimensional analyses, improving accuracy by incorporating factors such as the point spread function and non-axisymmetric components (Andredakis et al., 1995; Schade et al., 1995; Byun and Freeman, 1995).

The fitting of two-dimensional light profiles in galaxies has grown increasingly complex, requiring faster and more scalable algorithms. Two primary methods dominate: Markov Chain Monte Carlo (MCMC) and gradient-based approaches. MCMC-based programs, such as profit (Robotham et al., 2017) and autogalaxy (Nightingale et al., 2023), iteratively sample to approximate the target distribution, with fitted values and uncertainties derived from the Markov chain’s statistics. Conversely, gradient-based methods, utilized by tools like galfit (Peng et al., 2010) and imfit (Erwin, 2015), offer quicker solutions by directly navigating to the loss function’s local minima. However, this speed compromises the exploration of parameter space and the robustness of uncertainty estimates, a trade-off not present in MCMC-based approaches (Chen et al., 2016).

With the expanding volume of data from astronomical surveys such as the Sloan Digital Sky Survey (SDSS, York et al., 2000), Dark Energy Survey (DES, Collaboration: et al., 2016), Chinese Survey Space Telescope (CSST, Zhan, 2011), Legacy Survey of Space and Time (LSST, Abell et al., 2009), and EUCLID (Laureijs et al., 2011), significant challenges emerge in galaxy profile fitting tasks. Both gradient-based and MCMC approaches face difficulties in processing the vast number of galaxies efficiently.

To address these challenges, deep learning strategies are increasingly being employed. For instance, deeplegato (Tuccillo et al., 2018) and galnets (Li et al., 2022; Qiu et al., 2023) train convolutional neural networks to predict profile-fitting parameters from simulated galaxy images. These deep learning approaches establish an implicit link between galaxy images and their profile parameters, substantially accelerating the process of parameter estimation (see also gamornet, Ghosh, 2019). However, these methods cannot be easily adapted to new data domains due to the well-known ’black box’ problem (Castelvecchi, 2016), which means that the implicit mappings generated by the models usually lack clear interpretations.

Another method to address speed challenges involves parallel gradient-based methodology, which relies exclusively on GPU-accelerated matrix parallel computation. This is supported by popular frameworks such as pytorch (Paszke et al., 2019), tensorflow (Abadi et al., 2016), and jax (Bradbury et al., 2018), enabling the incorporation of automatic differentiation. Applications of this approach range from galaxy kinematic modeling (Bekiaris, 2017) to cosmological N-body simulations (Modi et al., 2021).

To speed up galaxy profile fitting, bringing it to the LSST-era and contribute to the literature on galaxy formation and evolution, we introduce galmoss. We aim to harness the advantages of both traditional and contemporary approaches by providing a framework that enables seamless parallelization while preserving the interpretability of more conventional methods. galmoss processes batches of galaxy photometric images as multidimensional arrays, facilitating efficient fitting on GPUs. This tensor-based approach permits dynamic management of data sizes during computations through the adjustment of batch sizes and the optimization of GPU utilization. For the fitting process, galmoss makes use of gradient descent due to its relatively low computational and memory requirements. Finally, galmoss provides two methods for uncertainty estimation. In the realm of galaxy profile fitting, astrophot (Stone et al., 2023) leverages GPU acceleration and automatic differentiation in its serial optimization process while preserving physical interpretability. While our package focuses on processing numerous images of individual galaxies, astrophot is optimized for fitting multiple galaxies within a single, crowded image.

This paper is structured as follows: Section 2 provides an overview of the galmoss workflow. Section 3 describes the data that we use to showcase the software’s performance. Section 4 details the generation of model images, using the Sérsic profile as an example, and introduces a set of built-in profiles along with advanced applications. Section 5 discusses the fitting process post-image generation. Finally, Section 6 showcases galmoss’ performance through case studies, with concluding remarks presented in Section 7.

2 General workflow of galmoss

Figure 1: The galmoss workflow: from parameter initialization and image processing to the fitting procedure and final data product generation, including the fitted model and residuals, all stored in FITS format.

Figure 1 shows the workflow of galmoss. It begins by reading a vector $p=\{p_G,p_L\}$ consisting of user-defined initial parameter values, along with the data image, sigma image (data uncertainty), mask image, and Point Spread Function (hereafter PSF) image. The first subset $p_G$ describes the geometric properties, which define the central position of the ellipse’s isophote, its major/minor axes, and the position angle. The second subset $p_L$ reads the initial values for the galaxy’s radial brightness profile $\zeta (𝐫)$.

With these initial values established, galmoss generates a series of image models following chosen profiles. These models are then convolved with the PSF of the image and augmented by adding a mean sky value to all pixels. This process results in the image model $M_{i,j}$, which is then compared to the image data $D_{i,j}$ and corresponding sigma images ${\sigma }_{i,j}$. The fitting process employs a ${\chi }^2$ loss function along with a mask ${\mathrm{\Omega }}_{i,j}$ to determine which pixels are included in the optimization. Internally, galmoss employs a gradient descent method embedded in pytorch.

After completing the fitting process, the package saves the data products, including the image block and fitted parameter values, to the designated output path. The image block files, saved in the FITS format, comprise the original image, the best-fit model image, and the residuals. With each component stored in a given FITS Head Data Unit (HDU), in the case of fitting multiple components (e.g., bulge + disc), each component can also be accessed individually alongside other data products.

3 Data

For showcasing galmoss capabilities, we gathered galaxy images ($128\times 128$ pixels) and PSF images ($41\times 41$ pixels) from the SDSS catalog. Our data selection was guided by the availability of independent galaxy parameters, as provided by the MaNGA PyMorph Photometric Value Added catalog (MPP-VAC-DR17, Domínguez Sánchez et al., 2022), encompassing 10,293 galaxies from the final MaNGA release (Bundy et al., 2014). The selected sample predominantly consists of galaxies with an average redshift z $\sim$ 0.03 and features a relatively uniform distribution of stellar mass (Collaboration, 2022). The catalog offers photometric parameters obtained through two-dimensional surface profile fitting. These parameters include results from single Sérsic and combined Sérsic + Exponential models, calculated using the pymorph pipeline (Vikram et al., 2010). pymorph integrates sextractor with galfit, the former package streamlines the estimation of initial parameters and the generation of stamps, and the latter aims to fit the galaxy profile.

In our selection process, we excluded galaxies located at the edges of plates where obtaining a complete $128\times 128$ pixel image was not possible. We also left out galaxies that failed in the single Sérsic fitting process within the MPP-VAC-DR17 catalog. This resulted in a final sample of 8,289 galaxies. We used only g-band images to ensure a good signal-to-noise ratio throughout the sample. Following pymorph, we generated masked images and initial values using sextractor (Bertin and Arnouts, 1996).

4 Image Generation and built-in Profiles

For a given set of galaxy or profile parameters, galmoss generates model images that represent the intensity or surface brightness of each pixel, based on the selected profile. galmoss provides six built-in radial profiles: Sérsic, Exponential disk, Ferrer, King, Gaussian, and Moffat, as well as a Flat sky model (with the option to include more complex sky models). These radial profiles are transformed into two-dimensional images by assuming circular symmetry. This process effectively projects the one-dimensional profile over a circular area in two dimensions to create the final image. To illustrate this projection process more concretely, we use the Sérsic profile as an example, given its widespread application in modeling galaxy profiles.

Furthermore, galmoss supports the combination and integration of different profiles. Detailed instructions on the practical use of the code can be found in the B, C and D.

4.1 Sérsic

Figure 2: g-band galmoss fitting result of J162123.19+322056.4, showing a good fitting quality. The left four images show the model, model residuals (data-model), data, and the residual distribution. On the right, a one-dimensional projection is shown, where the upper panel represents the model distribution with a black line and the orange filled markers with error bars represent the galaxy image’s flux distribution and its uncertainty. The lower panel represents the flux residual with blue filled markers with error bars.

The Sérsic profile (Sersic, 1968) is commonly used to fit the surface brightness distribution of elliptical galaxies, as well as the disk and bulge components of other galaxy types. The intensity $I(r)$ as a function of the radius $r$ is given by

$$I(r)=I_{\mathrm{e}}\exp \left\{-{\nu }_n\left[{\left(\frac{r}{r_{\mathrm{e}}}\right)}^{\frac{1}{n}}-1\right]\right\},$$ (1)

where $I_{\mathrm{e}}$ denotes the surface brightness at the effective radius $r_{\mathrm{e}}$, which encompasses half of the total profile flux. The Sérsic index $n$ dictates the profile’s curvature, and the parameter ${\nu }_n$ is computed numerically from the inverse cumulative distribution function of the gamma distribution. The versatility of the Sérsic profile lies in its ability to model a wide range of galaxy morphologies by adjusting the Sérsic index $n$. For example, setting $n=4$ transforms it into the de Vaucouleurs profile, known as the $r^{1/4}$ profile. An exponential disk profile is obtained with $n=1$, and $n=0.5$ results in a Gaussian profile, enabling users to explore these diverse profiles by simply varying $n$.

In a one-dimensional context, the radius $r$ refers to the distance from the galaxy center. In two dimensions, given a galaxy profile center $(x_{\mathrm{c}},y_{\mathrm{c}})$, the position angle of the ellipse profile $\theta$, and the axis ratio $q$, the circular radius $r$ is defined by (see e.g. Robotham et al., 2017):

$r$ $={\left[r_{\mathrm{maj}}^{(B+2)}+{\left({\displaystyle \frac{r_{\min }}{q}}\right)}^{(B+2)}\right]}^{\frac{1}{B+2}},$ (2)

where

$r_{\mathrm{maj}}$ $=\left|\cos (\theta )(x-x_{\mathrm{c}})\right|+\left|\sin (\theta )(y-y_{\mathrm{c}})\right|,$ (3)

$r_{\min }$ $=\left|-\sin (\theta )(x-x_{\mathrm{c}})\right|+\left|\cos (\theta )(y-y_{\mathrm{c}})\right|.$ (4)

The boxiness parameter $B$ introduces flexibility in the shape of the ellipse (Athanassoula et al., 1990), allowing for more general elliptical profiles. $B=0$ corresponds to a standard ellipse, while positive values result in more box-like shapes and negative values in more disc-like shapes. This flexibility is particularly useful in galaxy modeling, where a range of elliptical shapes can represent the diverse morphologies observed in galaxies.

Rather than directly using the surface brightness $I_{\mathrm{e}}$ in the Sérsic, magnitude $m$ and its corresponding zero point magnitude $m_0$ is used to specify the intensity level. This approach is consistent with methodologies employed by widely-used galaxy fitting tools such as galfit (Peng et al., 2010) and profit (Robotham et al., 2017), facilitating easier comparison with observational data. Both quantities are related as follow:

$$I_{\mathrm{e}}=\frac{r_{\mathrm{box}}{\nu }_n^{2n}{10}^{-0.4(m-m_0)}}{2\pi q{r}_{\mathrm{e}}^{2}n\mathrm{\Gamma }(2n)\exp ({\nu }_n)},$$ (5)

where

$$r_{\mathrm{box}}=\frac{\pi (2+B)}{2\beta (\frac{1}{2+B},\frac{1}{2+B})}.$$ (6)

In this formulation, $r_{\mathrm{box}}$ serves as a geometric correction factor in $I_{\mathrm{e}}$ to account for deviations from a perfect ellipse, influenced by the level of diskiness or boxiness. Specifically, when $B$ = 0, indicating a perfect ellipse, $r_{\text{box}}$ = 1, implying no geometric correction. The beta function $\beta (a,b)$, is calculated using the relationship $\beta (a,b)=\frac{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}{\mathrm{\Gamma }(a+b)}$, where $\mathrm{\Gamma }$ denotes the Gamma function.

Figure 2 presents the fitting results using the Sérsic profile for galaxy J162123.19+322056.4, which is one of the galaxies in our catalog. The quality of this model is evidenced by the residuals, highlighting the precision of the fit. For those interested in the implementation details, the code used for this analysis is provided in A.

4.2 Other available profiles

In this section, we briefly discuss the built-in profiles other than Sérsic. All these profiles but the Flat sky profile have circular radius $r$ that follows Equation 2 to Equation 4.

4.2.1 Exponential Disk Profile

In an exponential disk profile, the intensity is defined as follows:

$$I(r)=I_0\exp \left(-\frac{r}{r_{{}_{\mathrm{s}}}}\right),$$ (7)

where $I_0$ is the brightness of the profile’s surface in the center (at radial distance $r=0$), and $r_{{}_{\mathrm{s}}}$ is the disk scale-length.

The relation among the magnitude ($m,m_0)$, $r_{{}_{\mathrm{s}}}$ and $I_0$ is given by

$$I_0=\frac{r_{\mathrm{box}}{10}^{-0.4(m-m_0)}}{2\pi q{r}_{\mathrm{s}}^2},$$ (8)

where $q$ is the axis ratio, and $r_{\mathrm{box}}$ is as defined in equation 6.

4.2.2 Modified Ferrer Profile

The modified Ferrer profile, which is characterised by a nearly flat core and a rapidly truncated shape on the periphery, was originally proposed by Ferrers (1877) and later modified by Laurikainen et al. (2005). The intensity is defined as follows:

$$I(r)=I_0{\left[1-{\left(\frac{r}{r_{\mathrm{out}}}\right)}^{2-b}\right]}^a,$$ (9)

where $I_0$ is also the central surface brightness parameter, $r_{\mathrm{out}}$ is the outer truncation radius, and $b$ and $a$ are parameters governing the slopes of the truncation and core, respectively.

The relation between the magnitude ($m$, $m_0$) and the profile model parameters is given by

$$I_0=\frac{r_{\mathrm{box}}{10}^{-0.4(m-m_0)}}{\pi q{r}_{\mathrm{out}}^{2}a\beta (a,1+\frac{2}{2-b})},$$ (10)

where $q$ is also the axis ratio and $\beta$ is the beta function defined in equation 6. The modified Ferrer profile is particularly suitable to model the bar structure in galaxies (Blázquez-Calero et al., 2020; Dalla Bontà et al., 2018).

4.2.3 Empirical (Modified) King Profile

Since its initial illustration by King (1962), the empirical King profile has been extensively used to fit both galactic and extragalactic globular clusters (BAOlab222https://github.com/soerenslarsen/baolab presented by Larsen, 2014; Chies-Santos et al., 2007; Bonatto and Bica, 2008; Tripathi et al., 2023). The intensity of this profile is defined as follows:

(11)

Here, $I_0$ is also the central surface brightness parameter. The core radius, $r_{\mathrm{c}}$, signifies the scale at which the density starts to deviate from uniformity, while $r_{\mathrm{t}}$, the truncation radius, marks the boundary of the cluster. The global power-law factor, $\alpha$, dictates the rate at which the density declines with distance from the center. The concentration of stars in globular clusters can be defined using the parameters $r_{\mathrm{t}}$ and $r_{\mathrm{c}}$, denoted by $c=\log \left(\frac{r_{\mathrm{t}}}{r_{\mathrm{c}}}\right)$.

Figure 3: Panels comparing the measured Sérsic profile parameters for the selected $\sim$ 8,000 SDSS galaxies, along with $R^2$ estimation. The panels each display galaxy magnitude ($m$), position angle ($\theta$), axis ratio ($q$), effective radius ($r_e$) and Sérsic index ($n$). The transparency of the dots illustrates their concentration. The blue line shows where galfit result equals galmoss result. The histogram on top and right illustrates the distribution of measured values for galmoss and galfit, respectively. The comparison shows proper consistency in all panels.

4.2.4 Gaussian Profile

In a Gaussian profile, the intensity is defined as follows:

$$I(r)=I_0\exp \left(-\frac{r^2}{2{\sigma }^2}\right),$$ (12)

where $I_0$ is also the central surface brightness parameter, and $\sigma$ is the radial dispersion. In the galmoss implementation, the Full Width at Half Maximum (FWHM = 2.354$\sigma$) is used instead. The classical application of the Gaussian profile includes modelling a simple PSF and point sources.

4.2.5 Moffat Profile

The Moffat profile (Moffat, 1969), commonly used for modelling a realistic telescope PSF, defines its intensity as follows:

$$I(r)=I_0{\left[1+{\left(\frac{r}{r_{\mathrm{d}}}\right)}^2\right]}^{-n},$$ (13)

where

$$r_{\mathrm{d}}=\frac{\mathrm{FWHM}}{2\sqrt{2^{\frac{1}{n}}-1}}.$$ (14)

and $I_0$ also is the central density. The concentration index, $n$, dictates whether the distribution is more Lorentzian-like ($n$ =1) or Gaussian-like ($n\to \mathrm{\infty }$).

When considering the ellipse with axis ratio $q$ and boxiness parameter ($r_{\mathrm{box}}$ in equation 6), the relationship between the observed magnitude $m$ and profile model parameter is

$$I_0=\frac{r_{\mathrm{box}}(n-1){10}^{-0.4(m-m_0)}}{\pi q{r}_{\mathrm{d}}^2}.$$ (15)

Similar to the Gaussian profile, the Moffat profile can also be used to model point sources.

4.2.6 Flat Sky

Unlike other radial light profiles, in galmoss, the sky profile is controlled by the sky mean value $I_{\mathrm{sky}}$ across all pixels without a radial matrix:

$$I=I_{\mathrm{sky}}.$$ (16)

In galmoss, the only parameter needed in sky Profile is $I_{\mathrm{sky}}$. If a sky profile varies with position (radius), it can be included as a user-defined profile.

5 Image Fitting and Evaluation

Following the generation of model images, galmoss implements an extended $\chi 2$ likelihood and gradient descent optimization, leveraging python and pytorch functionalities. In addition, uncertainty estimation is available.

5.1 Maximum Likelihood Estimation

The galmoss package employs a ${\chi }^2$ likelihood:

$\ln ℒ=-{\displaystyle \frac{k}{2}}\ln 2-\ln \mathrm{\Gamma }(k/2)+\left({\displaystyle \frac{k}{2}}-1\right)\ln {\chi }^2-{\displaystyle \frac{{\chi }^2}{2}},$ (17)

${\chi }^2={\displaystyle \sum _{i=1}^{N_i}}{\displaystyle \sum _{j=1}^{N_j}}{\mathrm{\Omega }}_{i.j}{\displaystyle \frac{{(D_{i,j}-M_{i,j})}^2}{{\sigma }_{i,j}^2}},$ (18)

where $M_{i,j}$ represents the model at pixel $i,j$, and ${\sigma }_{i,j}$ represents the uncertainty of $D_{i,j}$. In the ${\chi }^2$ distribution, the degrees of freedom, $k$, is the difference between the number of pixels with galaxy flux that are not included in the mask and the number of model fit parameters.

Internally galmoss adopts gradient descent for parallel fitting. Though slower in convergence than other traditional methods (e.g. Levenberg-Marquardt) it has lower computational and memory demands, enabling efficient parallel fitting on GPUs.

5.2 Confidence Intervals

Galmoss employs two strategies to estimate uncertainties in galaxy images, primarily using a covariance matrix based on Gaussian distribution assumptions (Hogg et al., 2010). The parameter uncertainties are calculated from the covariance matrix’s diagonal, derived using the Jacobian matrix (e.g., Gavin, 2022) for computational efficiency. This method approximates uncertainties as ${\sigma }_{\mathrm{p}}=\sqrt{\mathrm{diag}{\left[J^{T}WJ\right]}^{-1}}$, where $J$ is the Jacobian matrix, and $W$ is diagonal with $W_{i,j}=1/{\sigma }_{i,j}^2$, offering a 1-$\sigma$ confidence level for the parameters. The second methodology for uncertainty estimation is bootstrapping, which involves resampling with replacement each galaxy images multiple times (typically 100 iterations) and refitting them. The uncertainty is given by ${\sigma }^2=\frac{1}{M}{\sum }_{i=1}^M{\left[m_i-\overline{m}\right]}^2$, where $\overline{m}$ is the mean value of the estimated parameter. Figure 4 illustrates a comparison of uncertainty estimation methods for a random subset of galaxies. We display the uncertainties determined through bootstrap analysis (red error bars) and those derived from the covariance matrix (blue error bars). Generally, bootstrap uncertainties encompass a broader range, suggesting a more conservative estimate.

Figure 4: A illustrative comparison of uncertainty estimations for fitted galaxy parameters using two methods: covariance matrix (depicted with blue error bars) and bootstrapping (depicted with red error bars), for a random selection of galaxies within our sample. Typically, the uncertainties derived from bootstrapping are observed to be larger than those obtained from the covariance matrix.

Figure 5: A speed-up performance evaluation between galmoss and galfit. Each line represents the speedup time normalized by the running time of galfit. The results show that galmoss’s speed surpasses galfit when the batch size slightly exceeds 10 and about 6 $\times$ faster when the batch size is 1,000. The benchmark was executed in a machine with the following specifications: CPU – 2.2GHz Intel Xeon Silver 4210; GPU – NVIDIA A100 (80GB); OS – Ubuntu Linux 18.04 64 bits; RAM – 3.9 TB.

6 Results

To evaluate the performance of galmoss against established methods, we performed a validation test by fitting single galaxy profiles from the same dataset and comparing the structural parameters, such as the Sérsic index and effective radius, with those obtained using galfit, as catalogued in the MPP-VAC-DR17. This comparison acts as a fundamental ’sanity check’ to verify the reliability of galmoss’s fitting results. Figure 3 shows overall consistency across key galaxy parameters—magnitude ($m$), position angle ($\theta$), axis ratio ($q$), effective radius ($r_{\mathrm{e}}$), and Sérsic index ($n$)—as quantified by the $R^2$ coefficient of determination.

The panels for magnitude, position angle, and axis ratio show a strong linear correlation, indicative of a high degree of alignment with galfit results. Specifically, in the position angle panel, objects that deviate significantly from the identity line are those fitted with large axis ratios ($q\sim 1$). This suggests that, in these cases, the light distributions are relatively insensitive to variations in the position angle ($\theta$).

The parameters $r_{\mathrm{e}}$ and $n$, known to be challenging to fit accurately (Trujillo et al., 2001), demonstrate a linear relationship, albeit with greater dispersion compared to other parameters. Notably, in the Sérsic index panel, galmoss values are generally lower than those derived from galfit at higher values of $n$. This discrepancy could be attributed to variations in image quality and the characteristics of the optimization algorithms used in pytorch. For an in-depth analysis of this observed bias, we direct the reader to Section E.

In addition to the aforementioned test, we conducted an independent speed comparison using galfit and galmoss using the same dataset and initial values. The fitting of 8,289 galaxies with galmoss completed in roughly 10 minutes—six times faster than with serial galfit. Figure 5 demonstrates this efficiency gain across different batch sizes, showing galmoss’ speed advantage becoming more pronounced as batch size increases, up to the limits of GPU capacity.

7 Conclusions

In this study, we introduce galmoss, an open-source software package specifically designed for fitting galaxy light profiles on large datasets, ideally suited for the LSST era. Built on the torch framework, galmoss provides efficient 2D surface brightness profile fitting for batches of galaxy images, with the added benefit of GPU acceleration. It features a user-friendly interface that allows for the easy definition of parameters and their ranges. Moreover, galmoss is capable of efficiently quantifying uncertainty through both covariance matrix and bootstrap methods. The package also supports the integration of new profiles, making it an adaptable and versatile tool for statistical analysis in galaxy structure studies.

Benchmark tests reveal that galmoss can achieve speedup gains of up to 6 times compared to galfit, depending on the sample size. One novel aspect of galmoss is its use of native parallelization to process multiple observational fields of single galaxies as elements of a single multidimensional tensor, thereby enhancing its speed through efficient PyTorch GPU-accelerated matrix calculations and memory usage (Paszke et al., 2019). This approach contrasts with astrophot, which focuses on larger fields with multiple galaxies or joint multi-band fitting. A current limitation of galmoss is the requirement to manually select profile models and initial parameters, along with its sub-optimal GPU utilization for small galaxy batches—areas for improvement in future versions.

galmoss is freely available at GitHub , Zenodo and listed in the Python Package Index . The readthedocs contains example usages, along with an overview of the package.

Appendix A How to use galmoss: single Sérsic case

Here, we demonstrate how to fit a single Sérsic profile to SDSS image data using the galmoss package.

First, we need to load the necessary packages.

⬇

1import Galmoss as gm

Next, we need to define the parameter objects and associate them with profile instances. The initial estimates of the galaxy parameters are provided by sextractor. Notably, we do not include the boxiness parameter in this simple example, despite its availability within the galmoss framework.

⬇

1# define parameter objects and profile

2

3sersic = gm.lp.Sersic(

4 cen_x=gm.p(65.43),

5 cen_y=gm.p(64.95),

6 pa=gm.p(-81.06, angle=True),

7 axis_r=gm.p(0.64),

8 eff_r=gm.p(7.58, pix_scale=0.396),

9 ser_n=gm.p(1.53, log=True),

10 mag=gm.p(17.68, M0=22.5)

11)

The comprehensive dataset object can be formulated utilising the image sets (galaxy image, sigma image, PSF image, mask image) together with the chosen profiles.

⬇

1dataset = gm.Dataset(

2 galaxy_index="J162123.19+322056.4",

3 image_path="./J162123.19+322056.4_image.fits",

4 sigma_path="./J162123.19+322056.4_sigma.fits",

5 psf_path="./J162123.19+322056.4_psf.fits",

6 mask_path="./J162123.19+322056.4_mask.fits"

7 mask_index=2,

8 img_block_path="./test_repo",

9 result_path="./test_repo"

10)

11

12dataset.define_profiles(sersic=sersic)

After initializing the hyperparameter during the fitting process, fitting could start. Subsequently, we run the uncertainty estimation process.

⬇

1fitting = gm.Fitting(dataset=dataset,

2 batch_size=1,

3 iteration=1000)

4fitting.fit()

5fitting.uncertainty(method="covar_mat")

When the fitting process is completed, the fitted results and the img_blocks are saved in corresponding path.

Appendix B How to use galmoss: bulge+disk case

Figure 6: g-band galmoss decomposition result of J100247.00+042559.8, showing a good fitting quality. The left six images show the total model, disk model, bulge model, data, model residuals (data-model), and the residual distribution. On the right, a one-dimensional projection is shown. The upper panel represents the model distribution with a black line, and the orange dots with error bars represent the galaxy image’s flux distribution and its uncertainty. The red dashed line and blue dotted line represent the distribution of the bulge component and the disk component, respectively. The lower panel represents the flux residual with blue dots with error bars.

Here, we demonstrate how to use a combination of two Sérsic profiles to make disk and bulge decomposition on SDSS image data using the galmoss package.

⬇

1import Galmoss as gm

Upon importing the package, the subsequent step entails defining parameter objects. To ensure that the center parameter within both profiles remains the same, it suffices to specify the center parameter once and subsequently incorporate it into various profiles.

⬇

1xcen = gm.p(65.97)

2ycen = gm.p(65.30)

For a quick start, we let the disk and bulge profile share the initial value from the sextractor, with an initial Sérsic index of 1 for the bulge component and 4 for the disk component.

⬇

1bulge = gm.lp.Sersic(

2 cen_x=xcen,

3 cen_y=ycen,

4 pa=gm.p(58.70, angle=True),

5 axis_r=gm.p(0.75),

6 eff_r=gm.p(4.09, pix_scale=0.396),

7 ser_n=gm.p(4),

8 mag=gm.p(17.97, M0=22.5)

9)

10

11disk = gm.lp.Sersic(

12 cen_x=xcen,

13 cen_y=ycen,

14 pa=gm.p(ini_value=58.70, angle=True),

15 axis_r=gm.p(0.75),

16 eff_r=gm.p(ini_value=4.09, pix_scale=0.396),

17 ser_n=gm.p(ini_value=1),

18 mag=gm.p(ini_value=17.97, M0=22.5)

19)

Compared to the single profile case, we only need to change the code snippet of profile definition. We choose to use bootstrap to calculate the uncertainty here.

⬇

1dataset = gm.Data_Box(

2 galaxy_index="J100247.00+042559.8"],

3 image_path="./J100247.00+042559.8_image.fits",

4 sigma_path="./J100247.00+042559.8_sigma.fits",

5 psf_path="./J100247.00+042559.8_psf.fits",

6 mask_path="./J100247.00+042559.8_mask.fits"

7 img_block_path="./test_repo",

8 result_path="./test_repo"

9)

10

11dataset.define_profiles(bulge=bulge, disk=disk)

12fitting = gm.Fitting(dataset=dataset,

13 batch_size=1,

14 iteration=1000)

15fitting.fit()

16fitting.uncertainty(method="bstrap")

Figure 6 show the decomposition result, with the residual demonstrating a high quality of fitting.

Appendix C Example code for definition of a set of galaxy profiles

Here, we provide an example that shows how to define a combination of profiles and extract the model image from the image function for an initial review. Figure 7 shows a model image featuring a central galaxy with a disk, bulge, and bar, each defined by two Sérsic profiles and a Ferrer profile, respectively. Additionally, there is a side galaxy defined by a Sérsic profile, two point sources defined by Moffat profiles, and an open cluster defined by a King profile. All of these are convolved with a PSF image produced by a Gaussian profile.

⬇

1import numpy as np

2import Galmoss as gm

3import torch_optimizer as optim

Upon importing the package, the first sub-model represents the central galaxy, comprising a disk, bulge, and bar.

⬇

1# The galaxy center

2xcen = gm.p(65.)

3ycen = gm.p(64.)

4

5# Define parameters and profiles

6Bisk = gm.lp.Sersic(

7 cen_x=xcen,

8 cen_y=ycen,

9 pa=gm.p(55., angle=True),

10 axis_r=gm.p(0.44),

11 eff_r=gm.p(13, pix_scale=0.396),

12 ser_n=gm.p(1),

13 mag=gm.p(17, M0=22.5),

14 box=gm.p(0.1)

15)

16

17Bulge = gm.lp.Sersic(

18 cen_x=xcen,

19 cen_y=ycen,

20 pa=gm.p(0., angle=True),

21 axis_r=gm.p(1),

22 eff_r=gm.p(6, pix_scale=0.396),

23 ser_n=gm.p(6),

24 mag=gm.p(21.5, M0=22.5)

25)

26

27Bar = gm.lp.King(

28 cen_x=xcen,

29 cen_y=ycen,

30 pa=gm.p(-40., angle=True),

31 axis_r=gm.p(0.2),

32 mag=gm.p(21.5, M0=22.5),

33 trunc_r=gm.p(5, pixScale=0.396),

34 trunc_a=gm.p(0.5),

35 trunc_b=gm.p(1.9),

36 box=gm.p(0.1)

37)

The second sub-model is a side galaxy.

⬇

1Sersic_s = gm.lp.Sersic(

2 cen_x=gm.p(45.43),

3 cen_y=gm.p(100.95),

4 pa=gm.p(-60., angle=True),

5 axis_r=gm.p(0.7),

6 eff_r=gm.p(13, pix_scale=0.396),

7 ser_n=gm.p(1.6),

8 mag=gm.p(17.5, M0=22.5),

9 box=gm.p(-0.6)

10)

The third sub-models are two point sources. They share the dispersion $\sigma$ = $1.515\mathrm{″}$ with the idealizsed PSF image, both of which are produced by a Gaussian profile.

⬇

1# Two point-souce

2P1 = gm.lp.Gaussian(

3 cen_x= gm.p(75),

4 cen_y= gm.p(80),

5 pa=gm.p(0., angle=True),

6 axis_r=gm.p(1),

7 inten=gm.p(0.1),

8 fwhm=gm.p(1.2,pix_scale=0.396),

9)

10

11P2 = gm.lp.Gaussian(

12 cen_x=gm.p(88.),

13 cen_y=gm.p(40.),

14 pa=gm.p(0., angle=True),

15 axis_r=gm.p(1),

16 inten=gm.p(0.1),

17 fwhm=gm.p(1.2, pix_scale=0.396),

18)

19

20# The idealized PSF

21PSF = gm.lp.Gaussian(

22 cen_x=gm.p(20),

23 cen_y=gm.p(20),

24 pa=gm.p(0., angle=True),

25 axis_r=gm.p(1),

26 inten=gm.p(1/(2*np.pi*((0.6/0.396)**2))),

27 fwhm=gm.p(1.2, pix_scale=0.396),

28)

The total model image without convolution is the sum of three sub-components (see Figure C.7(a)). To extract the model images, it is necessary to first define a grid, as grid generation occurs automatically only during the fit process within the fitting object.

⬇

1x = torch.linspace(0.5, 128-0.5, 128)

2y = torch.linspace(128-0.5, 0.5, 128)

3grid = torch.meshgrid(y, x, indexing=’ij’)

⬇

1bulge = Bulge.image_via_grid_from(grid, mode="initial_model")

2disk = Disk.image_via_grid_from(grid, mode="initial_model")

3bar = Bar.image_via_grid_from(grid, mode="initial_model")

4sersic_s = Sersic_s.image_via_grid_from(grid, mode="initial_model")

5p1 = P1.image_via_grid_from(grid, mode="initial_model")

6p2 = P2.image_via_grid_from(grid, mode="initial_model")

7result = bulge + disk + bar + sersic_s + p1 + p2

To mimic the effect of seeing, an idealized PSF image will be generated from a relatively small grid and then convolved with the model images (except for the point sources). The result is shown in Figure C.7(b).

⬇

1#Grid generation

2x_psf = torch.linspace(0.5, 40-0.5, 40)

3y_psf = torch.linspace(40-0.5, 0.5, 40)

4grid_psf = torch.meshgrid(x_psf, y_psf, indexing=’ij’)

5

6# PSF image

7psf = PSF.image_via_grid_from(grid_psf, mode="initial_model")

8

9#Convolving

10result_conv = torch.squeeze(torch.nn.functional.conv2d((bulge+disk+bar+sersic_s).reshape(1,1, 128,128), psf.reshape(1, 1, 40, 40), padding="same", groups=1)) + p1 + p2

(a) The Model image without convolving a PSF image.

(b) The Model image with a PSF image convolved.

Figure 7: The images of the model are provided for an initial review. The model consists of a central galaxy that has a disk, bulge, and bar, each defined by two Sérsic profiles and a Ferrer profile, respectively. Additionally, there is a side galaxy defined by a Sérsic profile, two point sources defined by Gaussian profile.

Appendix D User-defined profile

In Section 4.2.6, we introduced the sky profile. Here, we show how to integrate a more realistic sky profile as an example of a user-defined profile. In this new sky profile, the sky intensity is defined as follows:

$$I=I_0+k_{{}_{\mathrm{x}}}(x-x_0)+k_{{}_{\mathrm{y}}}(y-y_0),$$ (19)

where $(x_0,y_0)$ is the geometric center of the image, $I_0$ is the sky value at $(x_0,y_0)$, and $(k_{{}_{\mathrm{x}}},k_{{}_{\mathrm{y}}})$ describes the variation of the sky value along the x and y axes.

To integrate this new profile, we need to define a new class, which we have named $NewSky$. This class should include a function that guides the generation of the image model from the given parameters. Here are some caveats. Firstly, every profile class should inherit from the basic class LightProfile to access general light profile functions. Next, the parameters should be loaded into the profile class in the __init__ function, and set as attributions (e.g., self.sky_0 = sky_0).

⬇

1import galmoss as gm

2

3class NewSky(gm.LightProfile):

4 def __init__(self, sky_0, grad_x, grad_y):

5 super().__init__()

6 self.psf = False

7 self.sky_0 = sky_0

8 self.grad_x = grad_x

9 self.grad_y = grad_y

10

11 def image_via_grid_from(self,

12 grid,

13 mode="updating_model"):

14 return (self.sky_0.value(mode)

15 + self.grad_x.value(mode)

16 * (grid[1]-(grid[0].shape[1] + 1)/2)

17 + self.grad_y.value(mode)

18 * (grid[0]-(grid[0].shape[0] + 1)/2))

The equation for the profile is defined within the $image\mathrm{\_}via\mathrm{\_}grid\mathrm{\_}from$ function. Parameter values are extracted after the mode value, which has a default value of updating_model. This mode calls for values that are continuously updated throughout the fitting process and have already been broadcast to a suitable shape for multi-dimensional matrix calculations.

A newly defined profile can be used as follows:

⬇

1sky = NewSky(sky_0=gm.p(0.3),

2 grad_x=gm.p(2),

3 grad_y=gm.p(3))

Appendix E The fitted bias in high Sérsic index end

As shown in Figure 3, we observe a general bias in the high Sérsic index ($n$) end in the comparison results. Galaxies with large $n$ values (e.g., $n>6$) in the galfit results mostly exhibit lower values (e.g., ) in the galmoss results.

To further investigate the type of high $n$ galaxies, we use the morphological classification afforded by catalog MDLM-VAC-DR17. We identify 8,289 galaxies into 2,849 ETGs and 5,440 LTGs, and further identify ETGs into 1,891 Elliptical galaxies, 793 S0 galaxies and 165 undefined galaxies. We find that ETGs account for 34.37 $\%$ of the total galaxy samples and represent 88.78 $\%$ in the biased region. To investigate further, we plot a histogram for the $n$ value in E and S0 (Figure 8). In the plot, the histogram almost coincides between galmoss fitted results and galfit fitted results in S0 galaxies. However, the distribution of the histogram is very different between the galmoss fitted results and the galfit fitted results in elliptical galaxies. This result reveals that Elliptical galaxies contribute the most to the inconsistency in fitted $n$ values.

There are several reasons for the observed result. One of them is the limitation of observation. Most elliptical galaxies have a steep central core and a flat outer wing. Due to the poor signal-to-noise ratio in the flat wing, the ${\chi }^2$ is insensitive to the change of $n$, which acts as a small gradient during the optimization. Another possible reason is related to the properties of the optimizer. Most prevalent optimizers in pytorch use self-adaptive learning rates which reduce the learning rate under low gradient conditions to fasten convergence. As a result, the sensitivity of $n$ is further reduced at the high-value end compared to galfit, which adds to the difference in results between the two software in this case.

Figure 8: The histogram shows the distribution of $n$ values in E and S0 galaxies. Solid lines represent galmoss results, and dashed lines are for galfit. Purple and pink correspond to S0 and E, respectively. It is evident from the inconsistency between the values fitted by galmoss and galfit that Elliptical galaxies contribute most to the inconsistency in the high end of $n$ fitted results.