# Morphological Parameters and Associated Uncertainties for 8 Million Galaxies
in the Hyper Suprime-Cam Wide Survey

Aritra Ghosh (Airt/r oeghaSh) Department of Astronomy, Yale University, New
Haven, CT, USA  Yale Center for Astronomy and Astrophysics, New Haven, CT, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>C. Megan Urry Yale Center for
Astronomy and Astrophysics, New Haven, CT, USA  Department of Physics, Yale
University, New Haven, CT, USA  Aayush Mishra Department of Physics, Indian
Institute of Science Education and Research, Bhopal, India  Laurence
Perreault-Levasseur Department of Physics, Univesité de Montréal, Montréal,
Canada  Ciela Institute, Montréal, Canada  Mila - Quebec Artificial
Intelligence Institute, Montréal, Canada  Center for Computational
Astrophysics, Flatiron Institute, New York, NY, USA  Priyamvada Natarajan
Department of Astronomy, Yale University, New Haven, CT, USA  Department of
Physics, Yale University, New Haven, CT, USA  Black Hole Initiative, Harvard
University, Cambridge, MA, USA  David B. Sanders Institute for Astronomy,
University of Hawai‘i, Honolulu, HI, USA  Daisuke Nagai Department of
Astronomy, Yale University, New Haven, CT, USA  Yale Center for Astronomy and
Astrophysics, New Haven, CT, USA  Department of Physics, Yale University, New
Haven, CT, USA  Chuan Tian (田川) Yale Center for Astronomy and Astrophysics,
New Haven, CT, USA  Department of Physics, Yale University, New Haven, CT, USA
Nico Cappelluti Department of Physics, University of Miami, Coral Gables, FL,
USA  Jeyhan S. Kartaltepe Laboratory for Multiwavelength Astrophysics, School
of Physics and Astronomy,
Rochester Institute of Technology, Rochester, NY, USA  Meredith C. Powell
Kavli Institute of Particle Astrophysics and Cosmology, Stanford University,
Stanford, CA, USA  Amrit Rau  Department of Computer Science, Yale University,
New Haven, CT, USA  Ezequiel Treister Instituto de Astrofísica, Facultad de
Física, Pontificia Universidad Católica de Chile, Santiago, Chile

###### Abstract

We use the Galaxy Morphology Posterior Estimation Network (GaMPEN) to estimate
morphological parameters and associated uncertainties for $\sim 8$ million
galaxies in the Hyper Suprime-Cam (HSC) Wide survey with $z\leq 0.75$ and
$m\leq 23$. GaMPEN is a machine learning framework that estimates Bayesian
posteriors for a galaxy’s bulge-to-total light ratio ($L_{B}/L_{T}$),
effective radius ($R_{e}$), and flux ($F$). By first training on simulations
of galaxies and then applying transfer learning using real data, we trained
GaMPEN with $<1\%$ of our data-set. This two-step process will be critical for
applying machine learning algorithms to future large imaging surveys, such as
the Rubin-Legacy Survey of Space and Time (LSST), the Nancy Grace Roman Space
Telescope (NGRST), and Euclid. By comparing our results to those obtained
using light-profile fitting, we demonstrate that GaMPEN’s predicted posterior
distributions are well-calibrated ($\lesssim 5\%$ deviation) and accurate.
This represents a significant improvement over light profile fitting
algorithms which underestimate uncertainties by as much as $\sim 60\%$. For an
overlapping sub-sample, we also compare the derived morphological parameters
with values in two external catalogs and find that the results agree within
the limits of uncertainties predicted by GaMPEN. This step also permits us to
define an empirical relationship between the Sérsic index and $L_{B}/L_{T}$
that can be used to convert between these two parameters. The catalog
presented here represents a significant improvement in size ($\sim 10\times$),
depth ($\sim 4$ magnitudes), and uncertainty quantification over previous
state-of-the-art bulge+disk decomposition catalogs. With this work, we also
release GaMPEN’s source code and trained models, which can be adapted to other
data sets.

Extragalactic astronomy (506), Galaxies (573), Galaxy classification systems
(582), Astronomy data analysis (1858), Neural networks (1933), Convolutional
neural networks (1938)

††software: PyTorch (Paszke et al., 2019), Ignite (Fomin et al., 2020), MLFlow
(Chen et al., 2020), Numpy (Harris et al., 2020), Scipy (Virtanen et al.,
2020), Astropy (Robitaille et al., 2013; Price-Whelan et al., 2018), Pandas
(McKinney, 2010), Scikit-learn (Pedregosa et al., 2011), Matplotlib (Hunter,
2007),

##  1 Introduction

The morphology of galaxies has been shown to be related to various other
fundamental properties of galaxies and their environment, including galaxy
mass, star formation rate, stellar kinematics, merger history, cosmic
environment, and the influence of supermassive black holes (e.g., Bender et
al., 1992; Tremaine et al., 2002; Pozzetti et al., 2010; Wuyts et al., 2011;
Huertas-Company et al., 2016; Powell et al., 2017; Shimakawa et al., 2021;
Dimauro et al., 2022). Therefore, quantitative measures of the morphological
parameters for large samples of galaxies at different redshifts are of
fundamental importance in understanding the physics of galaxy formation and
evolution.

Distributions of morphological quantities alone can place powerful constraints
on possible galaxy formation scenarios. And when combined with other physical
quantities, they can provide key insights into evolutionary processes at play
or even reveal the role of new physical mechanisms that impact evolution
(e.g., Kauffmann et al., 2004; Weinmann et al., 2006; Schawinski et al., 2007;
van der Wel, 2008; Schawinski et al., 2014). However, such studies often
involve subtle correlations or hidden variables within strong correlations
that demand greater statistics and measurement precision than what has been
available in the preceding decades. An often-overlooked factor in such studies
has been the computation of robust uncertainties. The computation of full
Bayesian posteriors for different morphological parameters is crucial for
drawing scientific inferences that account for uncertainty and, thus are
indispensable in the derivation of robust scaling relations (e.g., Bernardi et
al., 2013; van der Wel et al., 2014) or tests of theoretical models using
morphology (e.g., Schawinski et al., 2014).

A quantitative description of galaxy morphology is typically expressed in
terms of its structural parameters — brightness, shape, and size — all of
which can be determined by fitting a single two dimensional analytic light
profile to the galaxy image (e.g., Van Der Wel et al., 2012; Tarsitano et al.,
2018). However, moving beyond single-component determinations by using
separate components to analyze galaxy sub-structure (e.g., disk, bulge, bar,
etc.) can provide us additional insights into the formation mechanisms of
these components: bulges, disks, and bars may be formed as a result of secular
evolution (e.g., Kormendy, 1979; Kormendy & Kennicutt, 2004; Genzel et al.,
2008; Sellwood, 2014) or due to the interaction of disk instabilities with
smooth and clumpy cold streams (e.g., Dekel et al., 2009a, b). In this sense,
contrary to what is often expected, bulges can also be formed without major
galaxy mergers.

Over the last decade, machine learning (ML) has been increasingly used for a
wide variety of tasks—from identifying exoplanets to studying black holes
(e.g., Hoyle, 2016; Kim & Brunner, 2017; Shallue & Vanderburg, 2018; Sharma et
al., 2020; Natarajan et al., 2021). Unsurprisingly, these algorithms have
become increasingly popular for determining galaxy morphology as well. (e.g.,
Dieleman et al., 2015; Huertas-Company et al., 2015; Tuccillo et al., 2018;
Ghosh et al., 2020; Hausen & Robertson, 2020; Walmsley et al., 2020; Cheng et
al., 2021; Vega-Ferrero et al., 2021; Tarsitano et al., 2022). The use of
these techniques has been driven by the fact that traditional methods of
analyzing morphologies—visual classification and template fitting—are not
scalable to the data volume expected from future surveys such as the Vera
Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al., 2019),
the Nancy Grace Roman Space Telescope (NGRST; Spergel et al., 2013), and
Euclid (Racca et al., 2016).

Most previous applications of ML to galaxy morphology produced broad,
qualitative classifications rather than numerical estimates of morphological
parameters. Tuccillo et al. (2018) did estimate parameters of single-component
Sérsic fits. However, they did not analyze galaxy sub-structures or provide
uncertainties. In order to address these challenges, in Ghosh et al. (2022),
we introduced The Galaxy Morphology Posterior Estimation Network (GaMPEN).
GaMPEN is a machine learning framework that estimates full Bayesian posteriors
for a galaxy’s bulge-to-total light ratio ($L_{B}/L_{T}$), effective radius
($R_{e}$), and flux ($F$). GaMPEN takes into account covariances between the
different parameters in its predictions and has been shown to produce
calibrated, accurate posteriors.

GaMPEN can also automatically crop input galaxy cutouts to an optimal size
before determining their morphology. This feature is critical given that
morphology-determination ML frameworks typically require input cutouts of a
fixed size, and thus, cutouts of a “typical” size often contain secondary
objects in the frame. By cropping out most secondary objects in the frame,
GaMPEN can make more accurate predictions over wide ranges of redshift and
magnitude.

Figure 1: The limiting absolute magnitudes probed by the Hyper Suprime-Cam
(HSC) Wide Survey and Sloan Digital Sky Survey (SDSS) at different redshifts.

In this paper, we use GaMPEN to estimate Bayesian posteriors for
$L_{B}/L_{T}$, $R_{e}$, and $F$ for $\sim 8$ million galaxies with $z\leq
0.75$ from the Hyper Suprime-Cam (HSC) Wide survey (Aihara et al., 2018). A
few recent works have studied the morphology of smaller subsets of HSC
galaxies: Shimakawa et al. (2021) classified $\sim 2\times 10^{5}$ massive HSC
galaxies into spiral and non-spiral galaxies, while Kawinwanichakij et al.
(2021) analyzed $\sim 1.5\times 10^{6}$ HSC galaxies using single-component
Sérsic fits. To date, the state-of-the-art morphological catalog that provided
bulge+disk decomposition parameters at low redshift has been that of Simard et
al. (2011), which used Sloan Digital Sky Survey (SDSS; York et al., 2000)
imaging to estimate morphological parameters of $\sim 1$ million $m<18$
galaxies, most with $z<0.2$. As Figure 1 shows, HSC imaging allows us to probe
much fainter magnitudes than SDSS with significantly better seeing (0.″85 for
HSC-W g compared to 1.″4 for SDSS g). The catalog presented in this paper
builds on Simard et al. (2011) by providing an order of magnitude increase in
sample size and probing four magnitudes deeper, to a higher redshift
threshold. Along with estimates of parameters, this catalog also estimates
robust uncertainties of the predicated parameters, which have typically been
absent from previous large morphological catalogs. This catalog represents a
significant step forward in our capability to quantify the shapes and sizes of
galaxies in our universe.

This paper also demonstrates that ML techniques can be used to study
morphology in new surveys, which do not have already-classified large training
sets available. Most previous works involving ML to study galaxy morphology
have depended on the availability of an extensive training set of real
galaxies with known properties from the same survey or a similar one. However,
if Convolutional Neural Networks (CNNs) are to replace traditional methods for
morphological analysis, we must be able to use them on new surveys that do not
have a morphological catalog readily available to be used as a training set.
In this paper, we demonstrate that by first training on simulations of
galaxies and then using real data for transfer learning (fine-tuning the
simulation-trained network), we can fully train GaMPEN, while having to label
$<1\%$ of our dataset. This work outlines an easy pathway to apply
morphological ML techniques to upcoming large imaging surveys.

In §2, we describe the HSC data used in this study, along with the simulated
two-dimensional light profiles. §3 and §4 provide a brief introduction to
GaMPEN and how we train it. In §5, we outline how we determined the
morphological parameters of $\sim 60,000$ HSC-Wide galaxies for transfer
learning and validating the results of our ML framework. §6 provides detailed
information on the accuracy of GaMPEN’s predictions. §7 compares our
predictions to that of other catalogs. We end with a summary and discussion of
our results in §8.

The full Bayesian posteriors for all $\sim 8$ million galaxies are being
released with the publication of this work. We are also releasing the source
code of GaMPEN, along with the trained models and extensive documentation and
tutorials on how to use GaMPEN. The public data release is described in
Appendix A.

##  2 Data

###  2.1 Hyper Suprime-Cam Data

We apply GaMPEN to g,r,i-band data from the Hyper Suprime-Cam (HSC) Subaru
Strategic Program Public Data Release 2 (PDR2; Aihara et al., 2019). The
Subaru Strategic Program, ongoing since 2014, uses the HSC prime-focus camera,
which provides extremely high sensitivity and resolving power due to the large
8.2 meter mirror of the Subaru Telescope.

Figure 2:  The filter used for each redshift bin is shown along with the
wavelength range sampled by each filter. The blue line shows where rest-frame
$450$ nm emission falls for redshifts labeled on the x-axis. As this figure
shows, the chosen filters allow us to consistently perform morphology
determination in the rest-frame g-band.

In order to have a large but uniform sample of HSC PDR2 galaxies, we focus on
the largest volume HSC survey, namely, the Wide layer, which covers $1400$
deg2 to the nominal survey depth in all filters and contains over 450 million
primary objects. To consistently perform morphology determination in the rest-
frame g-band across our entire sample, we use different filters for galaxies
at different redshifts, as shown in Figure 2. We use g-band for $z\leq 0.25$,
r-band for $0.25<z\leq 0.50$, and i-band for $0.50<z\leq 0.75$. Given HSC’s
typical seeing, all objects with sizes $\lesssim 5$ kpc cannot be resolved
beyond $z=0.75$. Therefore, given that the HSC-Wide light profiles of the
large majority of galaxies beyond this redshift will be dominated by the PSF,
we restrict this work till $z=0.75$ and will explore the application of GaMPEN
to higher redshift deeper HSC data in future work.

In order to select galaxies, we used the PDR2 galaxy catalog produced using
forced photometry on coadded images. The HSC data has been reduced using a
pipeline built on the prototype pipeline developed by the Rubin Observatory’s
Large Synoptic Survey Telescope’s Data Management system, and we refer the
interested reader to Bosch et al. (2018) for more details. We use the
extendedness_value flag only to select extended sources. The
extendedness_value flag relies on the difference between the Composite Model
(CModel) and PSF magnitudes to separate galaxies from stars, and contamination
(from stars) increases sharply for $m>23$ for median HSC seeing in the Wide
layer as outlined here111https://hsc-
release.mtk.nao.ac.jp/doc/index.php/stargalaxy-separation-2/. Thus, for each
redshift bin, we select galaxies with magnitude $m<23$ (in the appropriate
band). The full query to download the data in each redshift bin is available
in Appendix B.

Table 1: Data Characteristics  Sample | Redshift | Number | Imaging | Spectra
---|---|---|---|---
Name |  |  | Band |
Low-z | $z\leq 0.25$ | 820,566 | g | 7.4%
Mid-z | $0.25<z\leq 0.50$ | 2,937,871 | r | 2.4%
High-z | $0.50<z\leq 0.75$ | 4,247,208 | i | 1.6%
Figure 3: Redshift (top) and magnitude (bottom) distributions for the $\sim 8$
million galaxies used in this study. We used spectroscopic redshifts when
available and high-quality photometric redshifts otherwise. The spectroscopic
completeness of each sub-sample is shown in Table 1.

We use spectroscopic redshifts when available ($\sim 2.5\%$) and high-quality
photometric redshifts otherwise. The spectroscopic redshifts were collated by
the HSC Subaru Strategic Program (SSP) team from a wide collection of
spectroscopic surveys – zCOSMOS DR3 (Lilly et al., 2009), UDSz (Bradshaw et
al., 2013), 3D-HST (Momcheva et al., 2016), FMOS-COSMOS (Silverman et al.,
2015), VVDS (Le Fèvre et al., 2013), VIPERS PDR1 (Garilli et al., 2014), SDSS
DR12 (Alam et al., 2015), GAMA DR2 (Liske et al., 2015), WiggleZ DR1
(Drinkwater et al., 2010), DEEP2 DR4 (Newman et al., 2013), PRIMUS DR1 (Cool
et al., 2013), and VANDELS DR2 (Pentericci et al., 2018). The photometric
redshifts in HSC PDR2 were calibrated based on the above spectroscopic sample
and were calculated using the Direct Empirical Photometric code (DEmP; Hsieh &
Yee 2014), an empirical quadratic polynomial photometric redshift fitting
code, and Mizuki (Tanaka, 2015), a Bayesian template fitting code. For an
extended description, we refer the interested reader to Nishizawa et al.
(2020). To remove galaxies with unsecure photometric redshifts, we set the
photoz_risk_best parameter $<0.1$, which controls the risk of the photometric
redshift being outside the range $z_{\mathrm{true}}\pm
0.15(1+z_{\mathrm{true}})$, with 0 being extremely safe to 1 being extremely
risky. Using the cleanflags_any flag, we further excluded objects flagged to
have any significant imaging issues (in the relevant band) by the HSC
pipeline. This flag can be triggered by a wide range of issues; however, for
$\sim 80\%$ of cases, the flag was triggered by cosmic ray hits, as shown in
Appendix B. We checked and confirmed that none of the above cuts significantly
modified the redshift or magnitude distribution of our galaxy sample.

Figure 4: Four randomly chosen galaxy cutouts are shown here for each redshift
bin, with the object of interest at the center of each cutout. Note that most
of these cutouts have secondary objects in the frame, which can often cause ML
algorithms to produce spurious classifications. GaMPEN uses a Spatial
Transformer Network to crop most secondary objects out of the frame (see §3).

The above process resulted in the selection of $\sim 8$ million galaxies, with
$\sim 1$ million, $\sim 3$ million, $\sim 4$ million galaxies in the low-z,
mid-z, and high-z bins, as shown in Table 1. The magnitude and redshift
distribution of the data is shown in Figure 3. Using the HSC Image Cutout
Service222https://hsc-release.mtk.nao.ac.jp/das_cutout/pdr2/, we downloaded
cutouts for each galaxy with sizes of $30\arcsec$, $24\arcsec$, $16\arcsec$
for the low-, mid-, and high-z bins, respectively. These sizes are large
enough that they should capture all objects in the relevant bin. Using the
results of our light-profile fitting analysis on a sub-sample, as outlined in
§5, we expect $\sim 99.5\%$ of our sample to have a downloaded cutout size
$\geq 10\times$(size of the major axis of the galaxy). Figure 4 shows randomly
chosen examples of galaxy cutouts for each redshift bin.

###  2.2 Simulated Galaxies for Initial Training

As outlined later in §4.2, we perform GaMPEN’s initial round of training on
mock galaxy image cutouts, simulated to match HSC observations in the
appropriate band. To generate mock images, we used GalSim (Rowe et al., 2015),
the modular galaxy image simulation toolkit. GalSim has been extensively
tested and shown to yield very accurate rendered images of galaxies. We
simulated 150,000 galaxies in each redshift bin, with a mixture of both single
and double-component galaxies, in order to have a diverse training sample. To
be exact, $75\%$ of the simulated galaxies consisted of both bulge and disk
components, while the remaining $25\%$ had either a single disk or a bulge.

Table 2: Parameter Ranges of Simulated Galaxies

Sample | Component | Sérsic Index | Half-Light Radius | Flux | Axis Ratio | Position Angle
---|---|---|---|---|---|---
Name | Name |  | (arcsec) | (ADUs) |  | (degrees)
Low-z | Single-Component Galaxies
| 0.8 - 1.2 — 3.5 - 5.0aaThe single-component galaxies are equally divided between galaxies with a Sérsic index between 0.8 - 1.2 and galaxies with a Sérsic index between 3.5 - 5.0. | 0.1 - 5.0 | 30 - $1.35\times 10^{5}$ | 0.25 - 1.0 |  $-90.0$ \- $90.0$
|  |  | ($m_{g}\sim$ 14 - 23) |  |
Double-Component Galaxies
Disk | 0.8 - 1.2 | 0.1 - 5.0 | 0.0 - 1.0bbFractional fluxes are noted here. The bulge flux fraction is chosen such that for each simulated galaxy it is added with the disk flux fraction to give $1.0$. The total flux of the galaxies is varied between the values given in the top-row for each sample | 0.25 - 1.0 |  $-90.0$ \- 90.0
| Bulge | 3.5 - 5.0 | 0.1 - 3.0 | 1.0 - Disk. Comp.bbFractional fluxes are noted here. The bulge flux fraction is chosen such that for each simulated galaxy it is added with the disk flux fraction to give $1.0$. The total flux of the galaxies is varied between the values given in the top-row for each sample | 0.25 - 1.0 | Disk Comp. $\pm\,[0,15]$ccThe bulge position angle differs from the disk position angle by a randomly chosen value between $-15$ and $+15$ degrees.
Mid-z | Single-Component Galaxies
| 0.8 - 1.2 — 3.5 - 5.0aaThe single-component galaxies are equally divided between galaxies with a Sérsic index between 0.8 - 1.2 and galaxies with a Sérsic index between 3.5 - 5.0. | 0.1 - 3.0 | 30 - $3\times 10^{4}$ | 0.25 - 1.0 |  $-90.0$ \- $90.0$
|  |  | ($m_{r}\sim$ 15.5 - 23) |  |
Double-Component Galaxies
Disk | 0.8 - 1.2 | 0.1 - 3.0 | 0.0 - 1.0bbFractional fluxes are noted here. The bulge flux fraction is chosen such that for each simulated galaxy it is added with the disk flux fraction to give $1.0$. The total flux of the galaxies is varied between the values given in the top-row for each sample | 0.25 - 1.0 |  $-90.0$ \- 90.0
| Bulge | 3.5 - 5.0 | 0.1 - 2.0 | 1.0 - Disk. Comp.bbFractional fluxes are noted here. The bulge flux fraction is chosen such that for each simulated galaxy it is added with the disk flux fraction to give $1.0$. The total flux of the galaxies is varied between the values given in the top-row for each sample | 0.25 - 1.0 | Disk Comp. $\pm\,[0,15]$ccThe bulge position angle differs from the disk position angle by a randomly chosen value between $-15$ and $+15$ degrees.
High-z | Single-Component Galaxies
| 0.8 - 1.2 — 3.5 - 5.0aaThe single-component galaxies are equally divided between galaxies with a Sérsic index between 0.8 - 1.2 and galaxies with a Sérsic index between 3.5 - 5.0. | 0.1 - 2.0 | 30 - $2\times 10^{4}$ | 0.25 - 1.0 |  $-90.0$ \- $90.0$
|  |  | ($m_{i}\sim$ 16 - 23) |  |
Double-Component Galaxies
Disk | 0.8 - 1.2 | 0.1 - 2.0 | 0.0 - 1.0bbFractional fluxes are noted here. The bulge flux fraction is chosen such that for each simulated galaxy it is added with the disk flux fraction to give $1.0$. The total flux of the galaxies is varied between the values given in the top-row for each sample | 0.25 - 1.0 |  $-90.0$ \- 90.0
| Bulge | 3.5 - 5.0 | 0.1 - 1.5 | 1.0 - Disk. Comp.bbFractional fluxes are noted here. The bulge flux fraction is chosen such that for each simulated galaxy it is added with the disk flux fraction to give $1.0$. The total flux of the galaxies is varied between the values given in the top-row for each sample | 0.25 - 1.0 | Disk Comp. $\pm\,[0,15]$ccThe bulge position angle differs from the disk position angle by a randomly chosen value between $-15$ and $+15$ degrees.

Note. — The above table shows the ranges of the various Sérsic profile
parameters used to simulate mock HSC cutouts. $75\%$ of the simulated galaxies
have both disk and bulge components, and the remainder has either a disk or a
bulge component. All the simulation parameters are drawn from uniform
distributions.

For both the bulge and disk components, we used the Sérsic profile, and the
parameters required to generate the Sérsic profiles were drawn from uniform
distributions over ranges given in Table 2. For the disk and bulge components,
we allow the Sérsic index to vary between $0.8-1.2$ and $3.5-5.0$,
respectively. We chose to have varying Sérsic indices as opposed to fixed
values for each component in order to have a training set with diverse light
profiles. The parameter ranges for fluxes, and half-light radii are quite
expansive, such that the simulations represent most local galaxies (Binney &
Merrifield, 1998) at the appropriate redshift range. To obtain the
corresponding angular sizes for simulation, we placed the sample at the
appropriate redshift using Planck18 cosmology ($H_{0}=67.7$ km/s/Mpc, Planck
Collaboration et al., 2018) and the appropriate pixel scale.

Figure 5: Two stages of simulating an HSC galaxy. (Left): A randomly chosen
two-dimensional light profile generated by GalSim. (Right): The same image
after PSF convolution and noise addition. The white pixels represent (small)
negative values that arise from the process of noise addition.

To make the two-dimensional light profiles generated by GalSim realistic, we
convolved these with representative point-spread functions (PSFs) downloaded
from the HSC survey. We also added realistic noise using one-thousand
$2\arcsec\times 2\arcsec$ “sky objects” from the HSC PDR2 Wide field. Sky
objects are empty regions identified by the HSC pipeline that are outside
object footprints and are recommended for being used in blank-sky
measurements. For PSF convolution and noise addition, we follow the procedure
detailed in Ghosh et al. (2022) and refer the interested reader to that work
for more details. Figure 5 shows a randomly chosen simulated light profile and
the corresponding image cutout generated after PSF convolution and noise
addition. All the simulated galaxy images in each redshift bin were chosen to
have cutout sizes equal to their real data counterparts, as outlined in §2.1.

##  3 Brief Introduction to GaMPEN

The Galaxy Morphology Posterior Estimation Network (GaMPEN; Ghosh et al.,
2022) is a novel machine learning framework that can predict posterior
distributions for a galaxy’s bulge-to-total light ratio ($L_{B}/L_{T}$),
effective radius ($R_{e}$), and flux ($F$). In this section, we provide a
brief introduction to GaMPEN; however, for a complete understanding of
GaMPEN’s architecture and how it predicts posteriors, we refer the reader to
Appendix C and Ghosh et al. (2022).

The architecture of GaMPEN consists of an upstream Spatial Transformer Network
(STN) module followed by a downstream Convolutional Neural Network (CNN)
module. The upstream STN is used to automatically crop each input image frame
to an optimal size before morphology determination. Based on the input image,
the STN predicts the parameters of an affine transformation which is then
applied to the input image. The transformed image is then passed onto the
downstream CNN, which estimates the joint posterior distributions of the
morphological parameters. Because the transformation we use is differentiable,
the STN can be trained using standard backpropagation along with the
downstream CNN without any additional supervision.

The inclusion of the STN in the framework greatly reduces the amount of time
spent on data pre-processing as we do not have to worry about making cutouts
of the “proper size”—an inherently difficult task since most of the galaxies
in our sample have not been morphologically analyzed before. More importantly,
the STN greatly reduces the chance of spurious results by cropping most
secondary objects out of frame (see §6.2).

Figure 6: Diagram outlining the training (left) and posterior inference
(right) phases of the GaMPEN workflow. Training consists of feeding galaxies
(with pre-determined parameter values) through the STN and CNN modules,
minimizing the loss function using Stochastic Gradient Descent. During this
process, we re-scale the variables described in the text, and return them to
the original variable space during inference. After the STN+CNN networks are
trained, the posterior inference step consists of 500 forward passes with
dropout enabled for each galaxy image. We draw a sample from the predicted
multivariate Gaussian distribution during each forward pass, and the
collection of these samples gives us the predicted posterior distribution.

Two primary sources of error contribute to the uncertainties in the parameters
predicted by GaMPEN. The first arises from errors inherent to the input
imaging data (e.g., noise and PSF blurring), and this is commonly referred to
as aleatoric uncertainty. The second source of error comes from the
limitations of the model being used for prediction (e.g., the number of free
parameters in GaMPEN, the amount of training data, etc.); this is referred to
as epistemic uncertainty.

For every given input image, GaMPEN predicts a multivariate Gaussian
distribution $\mathcal{N}(\bm{\mu},\bm{\Sigma})$, with mean $\bm{\mu}$ and
covariance matrix $\bm{\Sigma}$. Although we would like to use GaMPEN to
predict aleatoric uncertainties, the covariance matrix, $\bm{\Sigma}$, is not
known a priori. Instead, we train GaMPEN to learn these values by minimizing
the negative log-likelihood of the output parameters for the training set. The
covariance matrix here represents the aleatoric uncertainties in GaMPEN’s
predictions.

In order to obtain epistemic uncertainties, we use the Monte-Carlo Dropout
technique (Srivastava et al., 2014), wherein during inference, each image is
passed through the trained networks multiple times. During each forward pass,
random neurons from the network are removed according to a Bernoulli
distribution, i.e., individual nodes are set to zero with a probability, $p$,
known as the dropout rate.

The entire procedure used to estimate posteriors is summarized in Figure 6.
Once GaMPEN has been trained, we feed each input image, $\bm{\hat{X}}_{n}$,
500 times into the trained model with dropout enabled. During each iteration,
we collect the predicted set of
$\left(\hat{\bm{\mu}}_{n,t},\bm{\hat{\Sigma}}_{n,t}\right)$ for the $t^{th}$
forward pass. Then, for each forward pass, we draw a sample
$\bm{\hat{Y}_{n,t}}$ from the multivariate normal distribution
$\mathcal{N}\left(\bm{\hat{\mu}}_{n,t},\bm{\hat{\Sigma}}_{n,t}\right)$. The
distribution generated by the collection of all 500 forward passes,
$\bm{\hat{Y}_{n}}$, represents the predicted posterior distribution for the
test image $\bm{\hat{X}}_{n}$. The different forward passes capture the
epistemic uncertainties, and each prediction in this sample also has its
associated aleatoric uncertainty represented by $\bm{\hat{\Sigma}_{n,t}}$.
Thus the above procedure allows us to incorporate both aleatoric and epistemic
uncertainties in GaMPEN’s predictions.

##  4 Training GaMPEN

Since most of the galaxies described in §2.1 have not been morphologically
analyzed before, we devise a method to train GaMPEN with only $0.5\%$ of our
sample of real galaxies. In order to achieve this, we first train GaMPEN using
the simulated galaxies described in §2.2. Thereafter, we apply “transfer
learning”, wherein we fine-tune the already trained network using a small
sample of real galaxies, analyzed using GALFIT (Peng et al., 2002).

Figure 7: Diagram outlining the different stages of training GaMPEN. We first
train GaMPEN using simulated light profiles described in §2.2. Thereafter, we
fine-tune the simulation-trained framework using $0.5\%$ of our real data
sample, for which we pre-determined the morphological parameters using light-
profile fitting, as described in §5. Finally, we process all the $\sim 8$
million galaxies in our dataset through the trained GaMPEN framework to obtain
estimates of their morphological parameters and associated uncertainties.

The process of training GaMPEN consists of the following steps, summarized in
Figure 7.

  * •

Simulating galaxies corresponding to the desired data set (here, HSC-Wide).

  * •

Initial training of GaMPEN on those simulated images.

  * •

Fine-tuning GaMPEN using a small amount of the real dataset (here, $\sim
0.5\%$) for which morphological parameters have been pre-determined. This
process is known as transfer learning.

  * •

Testing the performance of the trained network on a test set of real galaxies
(here, $\sim 0.25\%$) for which morphological parameters have been pre-
determined.

  * •

Processing the remainder of the real data (here, $\sim 99\%$) through the
trained framework.

We describe each of the above steps in detail below.

###  4.1 Data Transformations

To make the training process more robust against numerical instabilities, we
transform the input images and target variables following the steps outlined
in Ghosh et al. (2022).

Since reducing the dynamic range of pixel values has been found to be helpful
in neural network convergence (e.g., Zanisi et al., 2021; Walmsley et al.,
2021; Tanaka et al., 2022), we pass all images in our dataset through the
arsinh function. For all the target variables, we first apply the logit
transformation to $L_{B}/L_{T}$ and log transformations to $R_{e}$ and $F$:

$\bm{Y_{n}^{\prime}}=f^{\prime\prime}(\bm{Y}_{n}^{\prime\prime})=\left(\log\frac{L_{B}/L_{T}}{1-L_{B}/L_{T}},\log
R_{e},\log F\right),$ (1)

where $\bm{Y}_{n}^{\prime\prime}=[{L_{B}/L_{T},R_{e},F}]$ is the target set of
variables before transformation and $f^{\prime\prime}$ is how we will refer to
the transformation in Equation 1. Next, we apply the standard scalar
transformation to each parameter (calibrated on the training data), which
amounts to subtracting the mean value of each parameter and scaling its
variance to unity. These two transformations ensure that all three variables
have similar numerical ranges and prevent variables with larger numerical
values from making a disproportionate contribution to the loss function.

Post training, during inference, we apply the inverse of the standard scalar
function (with no re-tuning of the mean or variance), followed by the inverse
of the logit and log transformations, $f^{\prime\prime-1}$, as indicated in
Figure 6. Besides transforming the variables back to their original space,
these final transformations also ensure that the predicted values always
conform to physically meaningful ranges ($0\leq L_{B}/L_{T}\leq 1$; $R_{e}>0$;
$F>0$).

###  4.2 Initial Training of GaMPEN on simulated galaxies

The purpose of training GaMPEN initially on simulations is two-fold. Firstly,
it greatly reduces the amount of real data needed for training the framework.
Secondly, since simulated galaxies are the only situation where we have access
to the “ground-truth” morphological parameters, they provide the perfect
dataset to assess GaMPEN’s typical accuracy for the different output
parameters.

In Ghosh et al. (2022), we extensively tested and reported on GaMPEN’s
performance on simulated HSC $z\leq 0.25$ g-band galaxies. Here, we extend
that to include simulated r-band $0.25<z\leq 0.5$ and i-band $0.50<z\leq 0.75$
galaxies. Out of the 150,000 galaxies simulated in each z-bin, we use $75\%$
to train the framework and $15\%$ as a validation set to determine the optimum
value for various hyper-parameters (such as learning rate, batch size, etc.).
Thereafter, we use the remaining $15\%$, which the framework has never
encountered before, to evaluate the performance of the trained framework.

We train GaMPEN by minimizing its loss function using Stochastic Gradient
Descent. The different hyper-parameters that need to be tuned are: the
learning rate (the step size during gradient descent), momentum (acceleration
factor used for faster convergence), the strength of L2 regularization (the
degree to which larger weight values are penalized), and batch size (the
number of images processed before weights and biases is updated). To choose
these hyper-parameters, we trained GaMPEN with different sets of hyper-
parameters and chose the ones that resulted in the lowest value for the loss
function on the validation set. The final hyper-parameters for the trained
models are given in Appendix D.

Figure 8: Histograms of residuals for simulated galaxies (in the test set)
across all three redshift bins. We define the residuals as the difference
between the true value and the most probable value predicted by GaMPEN. The
dashed vertical line represents $x=0$, denoting cases with perfectly recovered
parameter values. The mean ($\mu$), median ($\tilde{\mu}$), and standard
deviation ($\sigma$) of each residual distribution are listed in each panel.

In order to test the robustness of our simulation-trained framework, we
compare the predictions made by GaMPEN on the test set to the true values
determined from the simulation parameters. The results are similar across all
redshift bins and closely follow what we determined previously in Ghosh et al.
(2022). The histograms of residuals for GaMPEN’s output parameters are shown
in Figure 8 across all three redshift bins. Note that to make all three
parameters dimensionless, we report the (Residual $R_{e}$)/$R_{e}$, instead of
simply Residual $R_{e}$. For each histogram, we also show the mean ($\mu$),
median ($\tilde{\mu}$), and standard deviation $(\sigma)$. The mean and median
help demonstrate that the distributions are all centered around zero, and the
standard deviation indicates the value of the “typical error” made by the
framework (i.e., $68.27\%$ of the time, GaMPEN’s prediction errors are less
than this value).

Note that for the simulated data, the typical error made by GaMPEN increases
with redshift. This is expected given that in our simulations, galaxies in the
higher redshift bins are preferentially smaller, fainter, and have lower
signal-to-noise ratios than their lower redshift counterparts—thus, these
galaxies are harder to analyze morphologically for any image processing
algorithm. However, for all the parameters across all redshift bins, the
GaMPEN error is typically always $\lesssim 15\%$ for the simulated sample.

###  4.3 Transfer Learning using Real Data

Transfer Learning as a data-science concept has been around since the late
1990s (e.g., Blum & Mitchell, 1998) and involves taking a network trained on a
particular dataset and optimized for a specific task, and re-tuning/fine-
tuning the weights and biases for a slightly different task or dataset. In
Ghosh et al. (2020), we introduced the concept of training on simulated light-
profiles and then transfer learning using real data for galaxy morphology
analysis. Here, we employ the same idea, taking the networks trained on the
simulated galaxies and then employing transfer learning using $\sim 15,000$
real HSC-Wide galaxies in each redshift bin.

Figure 9: Magnitude and redshift distributions for all galaxies in each
redshift bin, plotted along with the galaxies selected for transfer learning
(before and after applying quality cuts, as described in §5). Note that we
plot density on the y-axis, not the number of samples. The total number of
galaxies used for transfer learning is $\sim 0.5\%$ of all the galaxies in our
dataset. The relative density of some magnitude bins is higher than others
(e.g., $18<m\leq 20$ for low-z) because they span a smaller range while having
roughly the same number of galaxies as the other bins.

In order to select a sample of galaxies to use for transfer learning, we start
with the galaxies summarized in Table 1. Note that what matters in training a
CNN is not matching the observed distributions of the simulation parameters;
rather, it is spanning the full range of those parameters with high fidelity.
Having too many of a particular type—even if that is the reality in real
data—can result in lower accuracy for minority populations (e.g., Ghosh et
al., 2020). As seen in Figure 3, the sample for all three redshift bins is
heavily biased towards fainter galaxies. To ensure that GaMPEN gets to train
on enough bright galaxies, we split the mid- and high-z samples into two sub-
samples: $m\leq 21$ and $m>21$. For the low-z sample, since it has a more
substantial tail towards the lower m values (compared to the mid- and high-z
samples), we use three sub-samples—$m\leq 18$, $18<m\leq 20$, $m>20$. We
select 20,000 galaxies from each redshift bin, making sure to sample equally
across the magnitude bins mentioned above. Thereafter, we determine their
morphological parameters using light-profile fitting, as described later in
§5. The magnitude and redshift distributions of the selected galaxies are
shown in Figure 9. As seen from the figure, the transfer learning sample has
sufficiently large numbers of examples from all parts of the parameter space.
This empowers us to optimize GaMPEN for the full range of galaxy morphologies.

Out of the 20,000 galaxies selected in each redshift bin, we use $75\%$ for
fine-tuning the simulation trained GaMPEN models and another $5\%$ for
selecting the various hyper-parameters to be used during the transfer learning
process. The remaining $20\%$ is used to evaluate the performance of the fully
trained GaMPEN frameworks in §6.

In order to artificially augment the number of galaxies being used for
transfer learning, we apply random rotations and horizontal/vertical flips on
the galaxies earmarked for training. This takes the effective number of
samples used for fine-tuning in each bin from $\sim 15,000$ to $\sim 90,000$.
Using the values determined by light-profile fitting as the labels, we fine-
tune the three GaMPEN models trained on simulations. We choose the values of
the different hyper-parameters based on the loss computed on the validation
set, and the final chosen hyper-parameters for transfer learning are reported
in Appendix D.

####  4.3.1 Fine-Tuning the Dropout Rate

Figure 10: The calculated percentile coverage probabilities for different
dropout rates for the low-z bin. Note that the coverage probabilities have
been averaged over the three output variables. The coverage probabilities are
defined as the percentage of the total test examples where the value
determined using light profile fitting lies within a particular confidence
interval of the predicted distribution. A dropout rate of $4\times 10^{-4}$
leads to coverage probabilities very close to their corresponding confidence
levels. A similar process of tuning in the mid-z and high-z bins leads to an
optimal dropout rate of $2\times 10^{-4}$ in both of them.

Aside from the hyper-parameters mentioned in §4.2, there is one more
adjustable parameter in GaMPEN—the dropout rate, which directly affects the
calculation of the epistemic uncertainties outlined in §3. On average, higher
dropout rates lead networks to estimate higher epistemic uncertainties. To
determine the optimal value for the dropout rate, during transfer learning, we
trained variants of GaMPEN with dropout rates from $10^{-3}$ to $10^{-5}$, all
with the same optimized values of momentum, learning rate, and batch size
mentioned in Appendix D.

To compare these models, we calculate the percentile coverage probabilities
associated with each model, defined as the percentage of the total test
examples where the parameter value determined using light profile fitting lies
within a particular confidence interval of the predicted distribution. We
calculate the coverage probabilities associated with the $68.27\%$, $95.45\%$,
and $99.73\%$ central percentile confidence levels, corresponding to the
$1\sigma$, $2\sigma$, and $3\sigma$ confidence levels for a normal
distribution. For each distribution predicted by GaMPEN, we define the
$68.27\%$ confidence interval as the region on the x-axis of the distribution
that contains $68.27\%$ of the most probable values of the integrated
probability distribution. To estimate the probability distribution function
from the GaMPEN predictions (which are discrete), we use kernel density
estimation, which is a non-parametric technique to estimate the probability
density function of a random variable.

We calculate the $95.45\%$ and $99.73\%$ confidence intervals of the predicted
distributions in the same fashion. Finally, we calculate the percentage of
examples for which the GALFIT-ed parameter values lie within each of these
confidence intervals. An accurate and unbiased estimator should produce
coverage probabilities equal to the confidence interval for which it was
calculated (e.g., the coverage probability corresponding to the $68.27\%$
confidence interval should be $68.27\%$). For every redshift bin, we choose
the dropout rate for which the calculated coverage probabilities are the
closest to their corresponding confidence levels. This leads to a dropout rate
of $4\times 10^{-4}$ for the low-z bin, and $2\times 10^{-4}$ for the mid- and
high-z bins.

As an example, we show in Figure 10 the coverage probabilities (averaged
across the three output variables) for different dropout rates for the low-z
sample. As can be seen, higher values of the dropout rate lead to GaMPEN over-
predicting the epistemic uncertainties, resulting in too high coverage
probabilities. In contrast, extremely low values lead to GaMPEN under-
predicting the epistemic uncertainties. For a dropout rate of $4\times
10^{-4}$, the calculated coverage probabilities are very close to their
corresponding confidence levels, resulting in accurately calibrated
posteriors.

##  5 Galfitting Galaxies for Transfer Learning & Validation

Figure 11:  Steps used in our light-profile fitting pipeline to determine
morphological parameters, for two representative galaxies in each redshift
bin. From left to right we show the input image, the mask generated by Source
Extractor, the model generated by GALFIT, and the residuals. Note that since
we do not explicitly model any Fourier bending modes or coordinate rotations,
we expect features like spiral arms to show up in the residuals, as depicted
in the second row.  Figure 12: Morphological parameters determined for HSC-
imaged galaxies using our light-profile fitting pipeline, versus morphological
parameters for the same galaxies determined by Simard et al. (2011) based on
SDSS imaging. The black dashed diagonal corresponds to perfect agreement. The
vertical line in the middle panel shows the median SDSS g-band seeing.

In this section, we describe a semi-automated pipeline that we developed and
used to determine the morphological parameters of $\sim 60,000$ galaxies
($\sim 20,000$ in each z-bin), which are used for transfer learning and to
test the efficacy of the trained GaMPEN frameworks.

In order to estimate the parameters, we use GALFIT, which is a two-dimensional
fitting algorithm designed to extract structural components from galaxy
images. However, before running GALFIT, we run Source Extractor (Bertinl,
1996) on all the cutout frames in order to obtain segmentation maps for each
input image. We use these segmentation maps to mask all secondary objects
present in the cutout frame during light-profile fitting. We also use the
Source Extractor estimates of various morphological parameters as the initial
starting guesses during light profile fitting. Lastly, we use the Source
Extractor estimates to pick a cutout size for each galaxy, which we set to be
ten times the effective radius estimated by Source Extractor.

Using GALFIT, we fit each galaxy with two Sérsic components:- a disk-component
with Sérsic index, $n=1$, and a bulge component with $3.5<n<5.0$. For each
galaxy, we perform two consecutive rounds of fitting:- first with some
constraints placed on the values of the different parameters and second round
with almost no constraints. This is to help the parameters converge quickly
during the initial round while still allowing the full exploration of the
parameter-space in the subsequent round. In the initial round, we constrain
the radius of each galaxy to be between 0.5 - 90 pixels (0.85 - 15.12 arcsec)
and the magnitude to be between $\pm 7.5$ of the initial value guessed by
Source Extractor. We also constrain the difference between the magnitudes of
the two components to be between $-7.5$ and $+7.5$
($0.001<L_{B}/L_{T}<0.999$), and the relative separation between the centers
of the two light profiles to be not more than 30 pixels ($5.04$ arcsec). These
constraints are much more expansive than what we expect galaxy parameters to
be in this z-range and are similar to what has been used for many previous
studies (e.g., Van Der Wel et al., 2012; Tuccillo et al., 2018). During the
second round of fitting, we only use the constraints on the Sérsic indices and
the relative separation between the two components. After fitting the two
components, we calculate the $R_{e}$ of each galaxy as the radius that
contains $50\%$ of the total light in the combined disk + bulge model. Figure
11 shows the image cutouts, masks, fitted models, and residuals for two
typical galaxies chosen from each redshift bin.

After the two rounds of fitting, we excluded galaxies for which GALFIT failed
to converge ($\sim 3.2\%$, $\sim 3.8\%$, and $\sim 2.9\%$ for the low-, mid-,
and high-z bins, respectively). Thereafter, we visually inspected the fits of
a randomly selected sub-sample of $\sim 300$ galaxies from each redshift bin.
The process of visually inspecting $\sim 900$ galaxies helped us to identify
three failure modes of our fitting pipeline: i) for some galaxies, GALFIT
assigned an extremely small axis ratio to one of the components leading to an
unphysical lopsided bulge/disk; ii) for some galaxies, the centroids of the
two fitted components were too far away from each other; iii) some galaxies
had residuals that were too large. Examples of each of these failure modes are
shown in Appendix E.

Table 3: GALFIT Quality Cuts  Sample | Criteria | Excluded %
---|---|---
Low-z | (problematic_value_flags == True) | 26.04%
OR (Sep. Dist. $>2^{\prime\prime}$)
OR (max_iters_flag == True
| AND Red. $\chi^{2}>2.5$) |
Mid-z | (problematic_value_flags == True) | 32.10%
OR (Sep. Dist. $>1^{\prime\prime}$)
OR (max_iters_flag == True
| AND Red. $\chi^{2}>1.25$) |
High-z | (problematic_value_flags == True) | 34.13%
OR (Sep. Dist. $>0^{\prime\prime}.7$)
OR (max_iters_flag == True
| AND Red. $\chi^{2}>1.25$) |

In order to get rid of these problematic fits from our GALFIT-ed sample, we
use a mixture of GALFIT flags and calculated parameters. Specifically, as
summarized in Table 3) and detailed below:

  1. a)

problematic_value_flags333Refer to Bullet 8 at
https://users.obs.carnegiescience.edu/peng/work/galfit/TOP10.html flags
galaxies that have an extremely small axis-ratio ($<0.1$) or radius ($<0.5$
pixels), or any other parameters that caused issues with numerical
convergence.

  2. b)

max_iters_flag flags galaxies for which GALFIT quit after reaching the maximum
number of iterations (100).

  3. c)

The reduced $\chi^{2}$ of the fit is poor.

  4. d)

The distance between the centers of the two fitted components is too large.

Table 3 shows the criteria used to exclude fits in each redshift bin. The
variation in the thresholds with redshift is to account for the fact that
galaxies at higher redshift are preferentially smaller, fainter, and have
lower signal-to-noise ratios. Figure 9 shows that the exclusion criteria do
not selectively exclude more galaxies from certain regions of the
magnitude/redshift parameter space compared to others. In order to determine
appropriate thresholds for the various flags above, we balanced excluding too
many galaxies against ensuring only good fits are included in the final
transfer learning dataset. The choice of thresholds is arbitrary to a certain
extent. In order to empower users to retrain GaMPEN using different criteria
for their own scientific analysis, we are making public the entire catalog of
GALFIT-ed values as outlined in Appendix A.

To the best of our knowledge, there is no published large catalog of
bulge+disk decomposition of HSC galaxies against which we can compare the
results of our light-profile fitting pipeline. However, Simard et al. (2011)
performed bulge+disk decomposition using g and r-band imaging from the Sloan
Digital Sky Survey (SDSS). We should note that HSC-Wide differs from SDSS in a
multitude of ways, with the most significant differences being in median
seeing [HSC-Wide: $0.\arcsec 79$ compared to SDSS: $1.\arcsec 4$ in the
g-band] and pixel scale [HSC: 0.168 arcsecs/pixel compared to SDSS: 0.396
arcsecs/pixel]. Thus, we do not expect our analysis to yield the exact same
results as that of Simard et al. (2011). However, given that a significant
portion of our low-z sample overlaps with that of Simard et al. (2011), it is
still useful to compare our results to that of Simard et al. (2011), as the
overall trends should agree.

Using an angular cross-match diameter of $0.15$ arcsec, we cross-matched our
low-z GALFIT sample with that of Simard et al. (2011) to obtain a sample of
$\sim 6500$ galaxies. Note that although our HSC sample extends to $g<23$, the
cross-matched galaxies are mostly $g<20$ due to the shallower depth of the
SDSS data. In Figure 12, we compare the results of our light-profile fitting
results with that of Simard et al. (2011). The figure shows galaxies in
hexagonal bins of roughly equal size, with the number of galaxies in each bin
represented according to the colorbar on the right. Note that we are using a
logarithmic colorbar to explore the full distribution of galaxies, down to 1
galaxy/bin. Although there are outliers present for all three parameters, and
the scatter of the relationship depends on the parameter, our GALFIT-derived
parameters strongly correlate with that of Simard et al. (2011). According to
the Spearman’s rank correlation test (see Dodge, 2008, for more details),
there is a positive correlation for all three variables, and the null
hypothesis of non-correlation can be rejected at extremely high significance
($p<10^{-200}$). The correlation coefficients obtained for $L_{B}/L_{T}$,
$R_{e}$, and $F$ are 0.85, 0.95, and 0.98, respectively.

Note that the higher scatter in the $L_{B}/L_{T}$ relation is expected, given
that bulge+disk decomposition involves constraining two different imaging
components simultaneously and is thus more sensitive to algorithmic
differences and the differences in imaging quality between SDSS and HSC
mentioned above. These differences in imaging quality, data-reduction
pipeline, and filters (Kawanomoto et al., 2018) might also explain the slight
offsets in the $R_{e}$ and flux measurements. It is also important to note
that most of the differences in $R_{e}$ measurements are at effective radii
values smaller than the median SDSS g-band seeing, as seen in the middle panel
of Figure 12.

##  6 Evaluating GaMPEN’s performance

After both rounds of training are complete, we apply the final trained GaMPEN
models to all the $\sim 8$ million galaxies outlined in §2.1. In this section,
we evaluate GaMPEN’s performance to assess the reliability of its predictions.

###  6.1 Inspecting the Predicted Posteriors

Figure 13: Examples of predicted posterior distributions for two randomly
chosen galaxies from each redshift bin. The blue shaded histograms show the
predictions from GaMPEN, and the solid blue lines show the associated
probability distribution functions estimated by kernel density estimation.
These are used to calculate the confidence intervals shown in the figure with
pink, yellow, and green shading. The mode (solid red line) shows the most
probable value of each morphological parameter. As expected, in most cases,
the GALFIT-ed value (dashed black line) lies within the $68.27\%$ confidence
interval.

Using the procedure outlined in §3 and Figure 6, we use the trained GaMPEN
frameworks to obtain joint probability distributions of all the output
parameters for each of the $\sim 8$ million galaxies. Figure 13 shows the
marginalized posterior distributions for six randomly selected galaxies; along
with the input cutouts fed to GaMPEN. All the predicted distributions are
unimodal, smooth, and resemble Gaussian/skewed-Gaussian distributions. For
each predicted distribution, the figure also shows the parameter space regions
that contain $68.27\%$, $95.45\%$, and $99.73\%$ of the most probable values
of the integrated probability distribution. We use kernel density estimation
to estimate the probability distribution function (PDF; shown by a blue line
in the figure) from the predicted values. The mode of this PDF is what we
refer to as the “predicted value” henceforth. The figure also demonstrates how
GaMPEN predicts distributions of different widths based on the input frame
(e.g., the $L_{B}/L_{T}$ distribution in the second row is wider compared to
the first row), and we will explore this in more detail in §6.4.

By design, GaMPEN predicts only physically possible values. This is especially
apparent in the $L_{B}/L_{T}$ column of rows 1 and 4 of Figure 13. Note that
to achieve this, we do not artificially truncate these distributions. Instead,
we use data transformations, as outlined in §4.1. This ensures that the
predicted $L_{B}/L_{T}$ values are always between 0 and 1. Similarly, we also
ensure that the $R_{e}$ and $F$ values predicted by GaMPEN are positive
through appropriate transformations.

While performing quality checks on the predicted posteriors of all the $\sim
8$ million galaxies, we noticed that sometimes GaMPEN predicts $R_{e}$ and F
values outside the parameter range on which it was trained (e.g., $m>23$). It
has been shown that while machine learning frameworks are excellent at
interpolation, one should be extremely cautious while trying to extrapolate
too much beyond the training set (e.g., Quionero-Candela et al., 2009; Recht
et al., 2019; Taori et al., 2020). Thus, for each redshift bin, we exclude
galaxies with predicted $R_{e}$ and $F$ values which are outside the upper and
lower bounds of the training set by more than 0.5 arcsecs or 0.5 mags,
respectively. This led to $\sim 6\%$, $\sim 2.5\%$, and $\sim 1.2\%$ of the
data being excluded in the low-, mid-, and high-z bins.

###  6.2 Evaluating the STN performance

Figure 14: Examples of the transformation applied by the STN to two randomly
selected galaxies from each redshift bin. The top row shows the input galaxy
images, and the bottom row shows the corresponding output from the STN. The
numbers in the top-left yellow boxes help correspond the output images to the
input images. As can be seen, the STN learns to crop most secondary objects
present in the input frame.

As can be seen from Rows 3, 4, and 5 of Figure 13, GaMPEN can accurately
predict morphological parameters even when the primary galaxy of interest
occupies a small portion of the input cutout, and secondary objects are
present in the input frame. This is primarily enabled by the upstream STN in
GaMPEN which, during training, learns to apply an optimal amount of cropping
to each input image.

Figure 14 shows examples of the transformations applied by the STN to randomly
selected galaxies in all three redshift bins. As can be seen, the STN crops
out most secondary galaxies present in the cutouts and helps the downstream
CNN to focus on the galaxy of interest at the center.

Figure 15: (Left): Galaxies in the low-z bin with the lowest values of $s$
(i.e., the most aggressive crops) (Right): Galaxies in the low-z bin with the
highest values of $s$ (i.e., the least aggressive crops). The $s$ parameter
denotes the fraction of the input image that was retained in the STN output.
As can be seen, the STN correctly learns to apply the most aggressive crops to
small galaxies; and the least aggressive crops to large galaxies.

To further validate the performance of the STN, we measured the amount of
cropping applied by the STN for all galaxies in the low-z bin. We chose the
lowest redshift bin for this test as it has the largest range of galaxy radii
among the different redshift bins. After that, we sorted all the processed
images based on the amount of cropping applied to each input image. In Figure
15, we show example images from our dataset with extremely high and extremely
low values of applied crops. The $s$ parameter shown in the lower left of each
panel denotes what fraction of the input image was retained in the STN
output—higher values of $s$ denote that a more significant fraction of the
input image was retained in the output image produced by the STN (i.e.,
minimal cropping). Figure 15 demonstrates that (without us having to engineer
this specifically), the STN correctly learns to apply the most aggressive
crops to the smallest galaxies in our dataset, and the least aggressive crops
to the largest galaxies

Thus, GaMPEN’s STN learns to systematically crop out secondary galaxies in the
cutouts and focus on the galaxy of interest at the center of the cutout. At
the same time, the STN also correctly applies minimal cropping to the largest
galaxies, making sure the entirety of these galaxies remains in the frame.

###  6.3 Comparing GaMPEN predictions to GALFIT predictions

Figure 16: Percentile coverage probabilities achieved on the test set shown
separately for each redshift bin. The leftmost set of bars in each panel shows
the coverage probabilities when averaged over the three output parameters, and
the right three sets of bars show the coverage probabilities for each
parameter individually. The mean coverage probability never deviates by more
than $4.5\%$ from the claimed confidence interval, and when considered for
each parameter separately, the coverage probability never deviates by more
than $8.7\%$. This demonstrates that GaMPEN produces well-calibrated accurate
uncertainties.

Out of the $60,000$ galaxies analyzed using GALFIT in §5, we use $80\%$ as the
training and validation sets. We use the remaining $20\%$, which the trained
GaMPEN frameworks have never seen, to evaluate the accuracy of the predicted
parameters. We refer to this as the “test set” henceforth.

In Figure 16, we show the coverage probabilities achieved by GaMPEN on the
test set. Note that in § 4.3.1, we tuned the dropout rate using the validation
set, whereas the values in Figure 16 are calculated on the test set. In the
ideal situation, they would perfectly mirror the confidence levels; (e.g.,
$68.27\%$ of the time, the true value would lie within $68.27\%$ of the most
probable volume of the predicted distribution). Clearly, the coverage
probabilities achieved by GaMPEN are consistently close to the claimed
confidence levels, both when averaged over the three output parameters, as
well as for each parameter individually. The mean coverage probability never
deviates by more than $4.5\%$ from the claimed confidence interval, and when
considered for each parameter individually, the coverage probability never
deviates by more than $8.7\%$ from the corresponding confidence interval.
Additionally, we note that even for the case for which the coverage
probabilities are most discrepant ($68\%$ flux confidence interval for the
mid-z model), the uncertainties predicted by GaMPEN are in any case
overestimates (i.e., conservative). If GaMPEN were used in a scenario that
requires perfect alignment of coverage probabilities, users could employ
techniques such as importance sampling (Kloek & van Dijk, 1978) on the
distributions predicted by GaMPEN. We note here that incorporating the
covariances between the predicted parameters into our loss function was key to
achieving simultaneous calibration of all three output variables.

Figure 17: The most probable parameter values predicted by GaMPEN for all
galaxies in the test set plotted against the values determined using GALFIT.
Galaxies are plotted in hexagonal bins of roughly equal size, and the number
of galaxies in each bin is represented according to the logarithmic colorbar
to the right of each panel. The top, middle, and bottom rows show the results
for the low-, mid-, and high-z bins, respectively. The dashed black $y=x$ line
represents the line of equality. Across all three redshift bins, values
predicted by GaMPEN closely mirror the values obtained using light-profile
fitting.

Having shown above that the posteriors predicted by GaMPEN are well
calibrated, we now investigate how close the modes (most probable values) of
the predicted distributions are to the parameter values determined using
light-profile fitting. Figure 17 shows the most probable values predicted by
GaMPEN for the test set plotted against the values determined using GALFIT in
hexagonal bins of roughly equal size. The number of galaxies is represented
according to the colorbar on the right. Note that we use a logarithmic
colorbar to visualize even small clusters of galaxies in this plane, down to 1
galaxy/bin. Across all three redshift bins, a large majority of all the
galaxies are clustered around the line of equality, showing that the most
probable values of the distributions predicted by GaMPEN closely track the
values obtained using light-profile fitting. The scatter obtained for
$L_{B}/L_{T}$ is larger compared to the other two output parameters, and we
explore that in more detail later in this section.

Figure 18: Distributions of residuals for all galaxies in the test set;
specifically, the differences between the values predicted by GaMPEN and those
obtained via light-profile fitting. The top, middle, and bottom rows show the
results for the low-, mid-, and high-z bins, respectively. The boxes in the
top-left corner of each panel show the mean ($\mu$), median ($\tilde{\mu}$),
and standard deviation ($\sigma$) of each residual distribution. The $\sigma$
of each distribution identifies the typical disagreement for each parameter
(e.g., for the low-z bin, in $68.27\%$ cases, the predicted magnitude is
within $\pm 0.41$ of the value determined by light profile fitting). The
dashed red vertical line marks $x=0$.

In Figure 18, we show the residual distribution for GaMPEN’s output parameters
in all three redshift bins. We define the residual for each parameter as the
difference between the most probable value predicted by GaMPEN and the value
determined using light profile fitting. The box in the upper left corner gives
the mean ($\mu$), median ($\tilde{\mu}$), and standard deviation ($\sigma$) of
each residual distribution. All nine distributions are normally distributed
(verified using the Shapiro Wilk test), and have $\mu\sim\tilde{\mu}\sim 0$.
The $\sigma$ of each distribution also identifies the typical disagreement for
each parameter (e.g., for the low-z bin, in $68.27\%$ cases, the predicted
$L_{B}/L_{T}$ value is within $\pm 0.17$ of the value determined by light
profile fitting). The residual $R_{e}$, when converted to physical units,
correspond to typical disagreements of 0.32 arcsecs, 0.14 arcsecs, and 0.14
arcsecs for the low-, mid-, and high-z bins, respectively. The $L_{B}/L_{T}$
residuals are mostly constant across the three redshift bins, while the
disagreement between GaMPEN and GALFIT predictions for $R_{e}$ and $F$
decrease slightly as we go from the low to the higher redshift bins. This
could be driven by the fact that the HSC median seeing becomes better as we
move from the g-band to the i-band (g-band: $0.\arcsec 79$; r-band: $0.\arcsec
75$; i-band: $0.\arcsec 61$). Additionally, lower redshift galaxies have
preferentially more resolved structural features (e.g., spiral arms), which
are not accounted for in our disk + bulge GALFIT decomposition pipeline. This
could lead to a higher disagreement between the GALFIT and GaMPEN predictions.

Figure 19: Residuals of the output parameters (difference between GaMPEN and
GALFIT predictions) plotted against the values predicted by GaMPEN for all
galaxies in the low-z test set. To make the y-axis dimensionless for all three
parameters, we plot the fractional $R_{e}$ residuals instead of absolute
values. This figure allows us to assign quality labels to GaMPEN’s predictions
based on the output values (e.g., flagging regions of the parameter space with
high levels of disagreement, as shown by the line-shaded region in the top-
left panel). See § 6.3 for details. The equivalent figures for the mid- and
high-z bins are shown in Appendix F.

Although Figures 17 and 18 indicate the overall agreement between GaMPEN and
GALFIT, they do not reveal the dependence of this agreement on location in the
parameter space. This is critical information as this identifies regions of
parameter space where GaMPEN agrees especially well or badly with light-
profile fitting results, so that future users can flag the reliability of
predictions in these regions or use results only from certain regions of the
parameter space for specific scientific analyses. Figure 19 shows the
residuals for the three output parameters for the low-z bin plotted against
the values predicted by GaMPEN. Note that in order to make the y-axis
dimensionless for all three parameters, we plot the fractional $R_{e}$
residuals instead of absolute values. As in Figure 17, we have split the
parameter space into hexagonal bins and used a logarithmic color scale to
denote the number of galaxies in each bin. The trends between the residuals
and different parameters are very similar across all three redshift bins.
Thus, to keep the main text concise, we have shown the plot for the low-z bin
here and shown the same plot for the mid- and high-z bins in Appendix F.

For most of the panels, the large majority of galaxies are clustered uniformly
around the black dashed line, $y=0$, which denotes the ideal case of perfectly
recovered parameters (assuming the GALFIT parameters are correct). There are a
few notable features on the top row, which depicts the $L_{B}/L_{T}$
residuals. In the top left panel, the $L_{B}/L_{T}$ residuals are highest near
the limits of $L_{B}/L_{T}$. We noticed the same effect when testing GaMPEN
with simulated galaxies in Ghosh et al. (2022), and refer to this as the “edge
effect” For $L_{B}/L_{T}$ values near the edges (i.e., when the disk/bulge
component completely dominates over the other component), precisely
determining $L_{B}/L_{T}$ is challenging for GaMPEN — in fact, this is
difficult for any image analysis algorithm). In some of these cases, GaMPEN
assigns almost the entirety of the light to the dominant component, resulting
in the streaks seen at the edges of the figure. Poor structural parameter
determination for galaxies with $L_{B}/L_{T}<0.2$ and $L_{B}/L_{T}>0.8$ have
also been independently observed in other studies using different algorithms
(e.g., Euclid Collaboration et al., 2022; Häußler et al., 2022)

In order to mitigate this, GaMPEN users can choose to transform GaMPEN’s
quantitative predictions in the region $0.2\geq L_{B}/L_{T}\geq 0.8$
(demarcated by shaded lines in Figure 19) to qualitative values such as
“highly bulge-dominated” ($L_{B}/L_{T}\geq 0.8$) or “highly disk-dominated”
($L_{B}/L_{T}\leq 0.2$). We followed a similar procedure in Ghosh et al.
(2022) and found the net accuracy of these labels to be $\gtrapprox 95\%$.

The top-middle panel of Figure 19 also shows that $L_{B}/L_{T}$ residuals are
higher for galaxies with smaller $R_{e}$. In other words, GaMPEN and GALFIT
systematically disagree more for galaxies with smaller sizes—and this effect
becomes more pronounced as the sizes become comparable to the seeing of the
HSC-Wide Survey (g-band: $0.79$ arcsec).

###  6.4 Inspecting the Predicted Uncertainties

Figure 20: Uncertainties predicted by GaMPEN for each parameter plotted
against the predicted values for the low-z test set. The $\sigma$ for each
parameter is defined as the width of the $68.27\%$ confidence interval. Note
that we plot fractional uncertainties for the radius in order to make the
y-axis dimensionless for all three rows. The line-shaded region in the top-
left panel shows the region where we recommend transforming quantitative
$L_{B}/L_{T}$ predictions to qualitative labels (see §6.3 for details).

The primary advantage of a Bayesian ML framework like GaMPEN is its ability to
predict the full posterior distributions of the output parameters instead of
just point estimates. Thus, we would expect such a network to inherently
produce wider distributions (i.e., larger uncertainties) in regions of the
parameter space where residuals are higher.

Figure 20 shows the uncertainties for the three predicted parameters plotted
against the predicted values for the low-z test set. We define the uncertainty
predicted for each parameter as the width of the $68.27\%$ confidence interval
(i.e., the parameter interval that contains $68.27\%$ of the most probable
values of the predicted distribution; see Fig. 13). The y-axis of the middle
row has been normalized so that all three panels show dimensionless fractional
uncertainties. The distributions of uncertainties look very similar across all
redshifts; thus, we have shown the uncertainty distributions for the mid- and
high-z bins in Appendix F.

In Figure 19, we saw that GaMPEN’s $L_{B}/L_{T}$ residuals are higher for
lower values of $R_{e}$. Here, we see that GaMPEN accurately predicts higher
$L_{B}/L_{T}$ uncertainties for lower values of $R_{e}$. This compensatory
effect is what allows GaMPEN to achieve the calibrated coverage probabilities
shown in Figure 16.

In the right column, we see that for all three parameters, the uncertainty in
GaMPEN’s predictions increase for fainter galaxies. This is in line with what
we expect and had seen in Ghosh et al. (2022)—morphological parameters for
fainter galaxies are more difficult to constrain compared to brighter galaxies
and thus should have higher uncertainties.

The top left panel of Figure 20 shows that GaMPEN is reasonably certain of its
predicted bulge-to-total ratio across the full range of values but appears
slightly more certain when $L_{B}/L_{T}\leq 0.2$ or $L_{B}/L_{T}\geq 0.8$. We
had also seen the same effect in Ghosh et al. (2022) with simulated galaxies.
We found that the smaller uncertainties at the limits corresponded to the
single-component galaxies, while for the double-component galaxies, the edge
effect is less pronounced (see Figure 15 of Ghosh et al. (2022)). Here, we are
seeing the same effect— GaMPEN’s uncertainties at the edges are systematically
lower for galaxies that can be described completely by only a disk or bulge
component. These lower uncertainties also contribute to the higher residuals
near the edges of $L_{B}/L_{T}$ (as wider distributions would reduce the
number of galaxies with high residuals at the edges). However, as can be seen
from the top left panel of Figure 20, most of the galaxies with very low
uncertainties lie in the region $0.2\geq L_{B}/L_{T}\geq 0.8$, where we
recommend transforming the quantitative $L_{B}/L_{T}$ predictions to
qualitative labels, as outlined in §6.3.

The results shown in this section outline the primary advantage of using a
Bayesian framework like GaMPEN— even in situations where the network is not
perfectly accurate, it can predict the right level of precision, allowing its
predictions to be reliable and well-calibrated.

###  6.5 Comparing GaMPEN’s Uncertainty Estimates to Other Algorithms

Figure 21: Percentile coverage probabilities for the $68.27\%$ confidence
interval obtained by GaMPEN on our HSC sample compared to coverage
probabilities obtained by various light-profile fitting algorithms on
simulated Euclid data (from Euclid Collaboration et al., 2022). The rightmost
set of bars shows the values calculated on the entire dataset, while the other
sets display values calculated on sub-samples of galaxies with specific
magnitude ranges (AB mag, shown on the x-axis). Compared to light-profile
fitting tools, the uncertainties predicted by GaMPEN are better calibrated by
$\sim 15-25\%$ overall and by as much as $\sim 60\%$ for the brightest
galaxies.

As §6.3 and §6.4 demonstrate, GaMPEN predicts well-calibrated uncertainties on
our HSC data. Previous studies (e.g., Haussler et al., 2007) have shown that
analytical estimates of errors from traditional morphological analysis tools
like GALFIT or GIM2D (Simard et al., 2002) are smaller than the true
uncertainties by $\geq 70\%$ for most galaxies.

Recently, Euclid Collaboration et al. (2022) reported coverage probabilities
obtained by five different light-profile fitting tools—Galapagos-2 (Häußler et
al., 2022), Morfometryka (Ferrari et al., 2015), ProFit (Robotham et al.,
2017), and SourceXtractor++(Bertin et al., 2020)—on simulated Euclid data. The
Euclid sample consisted of $\sim 1.5$ million galaxies ranging from $I_{E}\sim
15$ to $I_{E}\sim 30$, simulated at $0.\arcsec 1$ /pixel. The simulations
included analytic Sérsic profiles with one and two components, as well as more
realistic galaxies generated with neural networks. As we did not have access
to the simulated Euclid dataset, we could not test GaMPEN’s performance on the
same data. Instead, we compared GaMPEN’s coverage probabilities for the HSC
data set to those reported for the Euclid simulations. Although the latter is
significantly different from our HSC sample, coverage probabilities reflect
the ability of the predicted uncertainty to capture the true uncertainty and
are not necessarily correlated with accuracy (which often varies across
different data sets). Moreover, GaMPEN’s predicted uncertainties can be tuned
for specific data sets, as shown in Figure 10, which should only improve the
GaMPEN outcome. Therefore, the results presented by Euclid Collaboration et
al. (2022) allow us to perform a preliminary comparison of GaMPEN’s
uncertainty prediction to that of other algorithms.

Figure 21 shows the $68.27\%$ coverage probabilities achieved by GaMPEN on the
HSC data compared to values for the four light profile fitting codes on the
simulated Euclid dataset (averaging over the different structural parameters).
When considering all galaxies, GaMPEN’s uncertainties are at least $\sim
15-25\%$ better calibrated than the other algorithms. The differences are much
larger for brighter galaxies, suggesting that the uncertainties predicted by
these algorithms depend primarily on the flux of the object. In contrast,
GaMPEN’s uncertainty predictions remain well-calibrated throughout and are
better by as much as $\sim 60\%$ for the brightest galaxies. The severe under-
prediction of uncertainties seems to hold true even for Bayesian codes like
ProFit. It is likely that GaMPEN’s robust implementation of aleatoric and
epistemic uncertainties, along with a carefully selected transfer learning set
spanning the entire magnitude range (see Fig. 9), allows it to predict well-
calibrated uncertainties across a wide range of magnitudes.

##  7 Comparing Our Predictions to Other Catalogs

Figure 22: GaMPEN predictions plotted against values estimated by Simard et
al. (2011) for a cross-matched sample of $\sim 20,000$ galaxies with $z<0.2$
and $m<19$. The density of points in each histogram is represented according
to the logarithmic colorbar on the right. The black dots show the median y
values in bins of equal width along the x-axis, with the error bars depicting
the average $68.27\%$ and $95.45\%$ confidence intervals predicted by GaMPEN
in that bin. The dash-dotted vertical line in the middle panel shows the
median SDSS g-band seeing. Figure 23: GaMPEN predictions plotted against
values estimated by Kawinwanichakij et al. (2021) for a cross-matched sample
of $\sim 200,000$ galaxies with $0.5<z\leq 0.75$ and $m\leq 23$. Similar to
Figure 22, we show the median y-values and associated error bars in bins of
equal width along the x-axis. The dash-dotted vertical line in the middle
panel shows the median HSC i-band seeing.

In §6.3, we compared GaMPEN’s predictions to the values determined using our
light-profile fitting pipeline. In order to further assess the reliability of
GaMPEN’s predictions, we now compare our predictions to two other
morphological catalogs.

In Figure 22, we compare GaMPEN’s predictions to the fits of Simard et al.
(2011). Using an angular cross-match diameter of $0.\arcsec 15$, we cross-
matched our entire sample to that of Simard et al. (2011) to obtain an
overlapping sample of $\sim 20,000$ galaxies. A large majority of these
galaxies have $m<19.5$ and $z<0.2$. We have included error bars in this figure
to show the typical width of predicted distributions for different regions of
the parameter space. We binned the x-axis of each parameter into bins of equal
width and plotted the median y-value in each bin as a point, with the error
bars showing the average $68.27\%$ and $95.45\%$ confidence intervals for all
galaxies in that bin.

Despite the differences in SDSS and HSC imaging quality (with regards to
pixel-scale and seeing), our results largely agree with those of Simard et al.
(2011) within the ranges of predicted uncertainties. Using Spearman’s rank
correlation test, we obtain correlation coefficients of 0.87, 0.96, and 0.98
for $L_{B}/L_{T}$, $R_{e}$, and $F$, respectively. It is also interesting to
note that the widths of the error bars account for almost the entire scatter
in the distribution of points, showing the robustness of GaMPEN’s uncertainty
estimates.

In the $L_{B}/L_{T}$ panel, there is a cluster of points near
$L_{B}/L_{T}=0.8$, and this is due to the edge-effect described in §6.3.
GaMPEN and Simard et al. (2011)’s $R_{e}$ predictions become discrepant (even
when accounting for uncertainty) for $R_{e}\sim 1\arcsec$. This is not
unexpected given that these $R_{e}$ values are almost $0.\arcsec 4$ smaller
than the median SDSS seeing, depicted by the dash-dotted vertical line in the
middle panel of Figure 22. The slight offset seen in the magnitude estimates
is due to the differences in SDSS and HSC imaging quality, data-reduction
pipeline, and filters (Kawanomoto et al., 2018), given that this discrepancy
disappears when we compare our results to another catalog that uses HSC
imaging (see below).

In Figure 23, we compare our results to that of Kawinwanichakij et al.
(2021)444Catalog obtained from John D. Silverman, private communication.,
wherein the authors fitted single Sérsic light-profiles to $1.5\times 10^{6}$
HSC i-band galaxies using LENSTRONOMY (Birrer & Amara, 2018), a multipurpose
open-source gravitational lens modeling Python package. Following Figure 22,
we also show mean error bars in different bins of the fitted parameter in this
figure. To compare results in the same band, we cross-matched our high-z
sample with that of Kawinwanichakij et al. (2021) to obtain a sample of $\sim
200,000$ galaxies with $0.50<z\leq 0.75$ and i $\leq 23$. Note that the
magnitude discrepancy present in Figure 22 now disappears. GaMPEN and
LENSTRONOMY radius measurements are also in agreement within the limit of
uncertainties, with the mean trend deviating slightly from the $y=x$ diagonal
at $R_{e}$ values lower than the HSC median seeing, depicted by the dash-
dotted vertical line in Figure 23. It is also important to note that while we
define $R_{e}$ as the radius that contains $50\%$ of the total light of the
combined bulge + disk profile, Kawinwanichakij et al. (2021) defines it as the
semi-major axis of the ellipse that contains half of the total flux of the
best-fitting Sérsic model. The correlation coefficients for $R_{e}$ and
magnitude are 0.87 and 0.96, respectively.

Since Kawinwanichakij et al. (2021) used single-component fits, we cannot
compare our $L_{B}/L_{T}$ predictions with their catalog. However, the left
panel of Figure 23 shows the correlation between fitted Sérsic index ($n$) and
measured bulge-to-total light ratio and can be used empirically to convert one
measure into another. We find that, in line with expectations, higher Sérsic
indices generally correspond to higher values of $L_{B}/L_{T}$. A large
majority of galaxies with $n\leq 1.5$ have most of their light in the disk
component. Galaxies with $n\geq 3$ have a large fraction of their light in the
bulge component, although this fraction may vary from values as low as $40\%$
to $95\%$. Note that these trends largely agree with what was reported by
Simmons & Urry (2008).

##  8 Conclusions & Discussion

In this paper, we used GaMPEN, a Bayesian machine learning framework, to
estimate morphological parameters ($L_{B}/L_{T}$, $R_{e}$, $F$) and associated
uncertainties for $\sim 8$ million galaxies in the HSC Wide survey with $z\leq
0.75$ and $m\leq 23$. Our catalog is one of the largest morphological catalogs
and is the first publicly available structural parameter catalog for HSC
galaxies. It provides an order of magnitude more galaxies compared to the
current state-of-the-art disk+bulge decomposition catalog of Simard et al.
(2011) while probing four magnitudes deeper and having a higher redshift
threshold. This represents an important step forward in our capability to
quantify the shapes and sizes of galaxies and uncertainties therein.

We also demonstrated that by first training on simulations of galaxies and
then utilizing transfer learning using real data, we are able to train GaMPEN
using $<1\%$ of our total dataset for training. This is an important
demonstration that ML frameworks can be used to measure galaxy properties in
new surveys, which do not have already-classified large training sets readily
available. Our implemented two-step process provides a new framework that can
be easily used for upcoming large imaging surveys like the Vera Rubin
Observatory Legacy Survey of Space and Time, Euclid, and the Nancy Grace Roman
Space Telescope.

We showed that GaMPEN’s STN is adept at automatically cropping input frames
and successfully removes secondary objects present in the frame for most input
cutouts. Note that the trained STN framework can be detached from the rest of
GaMPEN, and can be used as a pre-processing step in any image analysis
pipeline.

By comparing GaMPEN’s predictions to values obtained using light-profile
fitting, we demonstrated that the posteriors predicted by GaMPEN are well-
calibrated and GaMPEN can accurately recover morphological parameters. We note
that the full computation of Bayesian posteriors in GaMPEN represents a
significant improvement over estimates of errors from traditional
morphological analysis tools like GALFIT, Galapagos, or ProFit. We
demonstrated that while the uncertainties predicted by GaMPEN are well
calibrated with $\lesssim 5\%$ deviation across a wide range of magnitudes,
traditional light-profile fitting algorithms underestimate uncertainties by
$\sim 15-60\%$ depending on the flux of the galaxy being analyzed. These well-
calibrated uncertainties will allow us to use GaMPEN for the derivation of
robust scaling relations (e.g., Bernardi et al., 2013; van der Wel et al.,
2014) as well as for tests of theoretical models using morphology (e.g.,
Schawinski et al., 2014).

GaMPEN’s residuals increase for smaller galaxies, but GaMPEN correctly
accounts for that by predicting correspondingly higher uncertainties for these
galaxies. GaMPEN’s $L_{B}/L_{T}$ residuals are also high when the bulge or
disk component completely dominates over the other component. GaMPEN’s
quantitative $L_{B}/L_{T}$ predictions for these galaxies ($0.2\geq
L_{B}/L_{T}\geq 0.8$) can be transformed into highly accurate qualitative
labels with $\gtrsim 95\%$ accuracy. We did not detect a decline in GaMPEN’s
performance for fainter galaxies, and GaMPEN’s residual patterns were also
fairly consistent across all three redshift bins used in this study.

In order to assess the reliability of our catalog using a completely
independent analysis, we compared GaMPEN’s predictions to those of Simard et
al. (2011) and Kawinwanichakij et al. (2021). Within the limit of
uncertainties predicted by GaMPEN, our results agree well with these two
catalogs. We noticed a slight discrepancy in the flux values determined using
SDSS and HSC imaging, which can be attributed to the differences in imaging
quality, data-reduction pipeline, and filters between the two surveys.
Comparing our $L_{B}/L_{T}$ predictions to Kawinwanichakij et al. (2021)’s
measured Sérsic indices also allowed us to study the correlation between these
two parameters, and this can also be used empirically to switch between these
two parameters. We found that although galaxies with $n\geq 3$ have a large
fraction of their light in the bulge component, this fraction can be anywhere
between $40\%$ to $95\%$.

Similar to other previous structural parameter catalogs (e.g., Simard et al.,
2011; Tarsitano et al., 2018), we did not explicitly exclude merging galaxies
from this catalog. We are currently working to incorporate a prediction flag
within GaMPEN that will flag merging/irregular galaxies and highly blended
sources so that they can be analyzed separately — we will demonstrate this in
future publications. However, we would like to note that for the redshift
ranges and magnitudes considered in this study, blending and mergers only
constitute a limited fraction of the total sample. Using the
m_(g|r|i)_blendedness_flag and m_(g|r|i)_blendedness_abs_flux parameters
available in HSC PDR2, we estimate that $\sim 10\%$, $\sim 4\%$, and $\sim
3\%$ of galaxies in the low-z, mid-z, and high-z redshift bins, respectively,
have nearby sources within the parent object footprint that can affect
structural parameter measurements. For an extended description of these flags,
we refer the interested reader to Bosch et al. (2018) and note that users can
choose to ignore/exclude these galaxies from their analysis by setting
different values/thresholds for these two parameters.

With this work, we are publicly releasing: a) GaMPEN’s source code and trained
models, along with documentation and tutorials; b) A catalog of morphological
parameters for the $\sim 8$ million galaxies in our sample along with robust
estimates of uncertainties; c) The full posterior distributions for all $\sim
8$ million galaxies. All elements of the public data release are summarized in
Appendix A.

Although we used GaMPEN here to predict $L_{B}/L_{T}$, $R_{e}$, and $F$, it
can be used to predict any morphological parameter when trained appropriately.
Additionally, although GaMPEN was used here on single-band images, we have
tested that both the STN and CNN modules in GaMPEN can handle an arbitrary
number of channels, with each channel being a different band. We defer a
detailed evaluation of GaMPEN’s performance on multi-band images for future
work. Additionally, GaMPEN can also be used to morphologically analyze
galaxies from other ground and space-based observatories. However, in order to
apply GaMPEN on these data sets, one would need to perform appropriate
transfer learning using data from the target dataset.

Finally, to give readers an estimate of GaMPEN’s run-time, we note that once
trained, it takes GaMPEN $\sim 1$ millisecond on a GPU and $\sim 150$
milliseconds on a CPU to perform a single forward pass on an input galaxy.
These numbers, of course, change slightly based on the specifics of the
hardware being used. However, as these numbers show, even with access to $\sim
1000$ CPUs or $\sim 10$ GPUs, GaMPEN can estimate full Bayesian posteriors for
millions of galaxies in just a few days. Therefore, GaMPEN is fully ready for
the large samples expected soon from Rubin-LSST and Euclid. Determinations of
structural parameters along with robust uncertainties for these samples will
allow us to characterize both morphologies as well as other relevant
properties traced by morphology (e.g., merger history) as a function of cosmic
time, mass, and the environment with unmatched statistical significance.

## acknowledgments

This material is based upon work supported by the National Science Foundation
under Grant No. 1715512.

CMU and AG would like to acknowledge support from the National Aeronautics and
Space Administration via ADAP Grant 80NSSC18K0418.

AG would like to acknowledge the support received from the Yale Graduate
School of Arts & Sciences through the Dean’s Emerging Scholars Research Award.

AG would like to acknowledge computing grants received through the Amazon
Cloud Credits for Research Program and the Yale Center for Research Computing
(YCRC) Research Credits Program. AG would also like to thank the Yale Center
for Research Computing and Yale Information Technology Services staff members
and scientists, especially Robert Bjorson and Craig Henry, for their guidance
and assistance in the vast amount of computation required for this project
that was performed on Yale’s Grace computing cluster and the Yale Astronomy
compute nodes.

We would like to thank John D. Silverman and Lalitwadee Kawinwanichakij for
making the Kawinwanichakij et al. (2021) catalog available to us. AG would
like to thank Tim Miller and Imad Pasha for helpful discussions.

PN gratefully acknowledges support at the Black Hole Initiative (BHI) at
Harvard as an external PI with grants from the Gordon and Betty Moore
Foundation and the John Templeton Foundation.

ET acknowledges support from ANID through Millennium Science Initiative
Program - NCN19_058, CATA-BASAL ACE210002 and FB210003, and FONDECYT Regular
1190818 and 1200495.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical
communities of Japan and Taiwan, and Princeton University. The HSC
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.
Funding was contributed by the FIRST program from Japanese Cabinet Office, 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.

This paper makes use of software developed for the Large Synoptic Survey
Telescope. We thank the LSST Project for making their code available as free
software at http://dm.lsst.org.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of
the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project
Office, the Max-Planck Society and its participating institutes, the Max
Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for
Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham
University, the University of Edinburgh, Queen’s University Belfast, the
Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory
Global Telescope Network Incorporated, the National Central University of
Taiwan, the Space Telescope Science Institute, the National Aeronautics and
Space Administration under Grant No. NNX08AR22G issued through the Planetary
Science Division of the NASA Science Mission Directorate, the National Science
Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos
Lorand University (ELTE) and the Los Alamos National Laboratory.

Based, in part, on data collected at the Subaru Telescope and retrieved from
the HSC data archive system, which is operated by Subaru Telescope and
Astronomy Data Center at National Astronomical Observatory of Japan.

##  Appendix A Data Access

GaMPEN’s source code is being publicly released along with the publication of
this article. Along with the source code, we are also releasing documentation,
tutorials, the trained models, the entire catalog of morphological predictions
as well as the estimated PDFs for different morphological parameters of the
$\sim 8$ million galaxies in our sample.

Links to the various components of the public data release mentioned above can
be accessed at:-

  * •

http://gampen.ghosharitra.com/

  * •

http://www.astro.yale.edu/aghosh/gampen.html

  * •

https://github.com/aritraghsh09/GaMPEN

  * •

https://gampen.readthedocs.io/en/latest/

We caution users to carefully read the Public Data Release Handbook available
at https://gampen.readthedocs.io/en/latest/Public_data.html to understand
various aspects of the data release before using the morphological catalogs
produced as a part of this paper.

Along with the source code and trained models, this public data release also
includes the various parameters obtained using the light profile fitting
pipeline described in §5.

##  Appendix B Additional Details About Data

As outlined in §2, we used the cleanflags_any parameter available as part of
HSC PDR2 to exclude objects flagged to have any significant imaging issues by
the HSC pipeline. The various triggers which contribute to the above flag, as
well as their prevalence among the galaxies which we excluded from the
analysis, are shown in Figure 24. As can be seen, $\sim 80\%$ of the triggers
are caused by cosmic ray hits/interpolated pixels.

In addition, the full SQL queries used to download the low-, mid-, and high-z
data are shown in Listings 25, 26, and 27. Note that after downloading the
data using these queries, we further excluded data based on the flags referred
to in the previous paragraph, as well as the quality of photometric redshift
estimates. For an extended description, please refer to §2. As noted in §8,
users may additionally choose to use the various blendedness flags available
in HSC PDR2 to further exclude merging/blended galaxies.

Figure 24: We exclude $\sim 20\%$ of our downloaded galaxies due to different
flags being triggered. The distribution of various flags that contribute to
galaxies being excluded from our sample is shown in this figure. As can be
seen, the large majority of exclusions are due to cosmic-ray hits (and hence,
interpolated pixels).

⬇

SELECT

object_id

, ra

, dec

FROM

pdr2_wide.smallcat

LEFT JOIN pdr2_wide.specz USING (object_id)

LEFT JOIN pdr2_wide.photoz_mizuki USING (object_id)

WHERE

((photoz_best > -1 AND photoz_best <= 0.25) OR (specz_redshift > -1 AND
specz_redshift <= 0.25))

AND (g_kronflux_mag < 23 OR g_cmodel_mag < 23)

AND g_extendedness_value = 1

;

Figure 25: Low-z Sample SQL query

⬇

\parSELECT

object_id

, ra

, dec

FROM

pdr2_wide.smallcat

LEFT JOIN pdr2_wide.specz USING (object_id)

LEFT JOIN pdr2_wide.photoz_mizuki USING (object_id)

WHERE

((photoz_best > 0.25 AND photoz_best <= 0.50) OR (specz_redshift > 0.25 AND
specz_redshift <= 0.50))

AND (r_kronflux_mag < 23 OR r_cmodel_mag < 23)

AND r_extendedness_value = 1

;

Figure 26: Mid-z Sample SQL query

⬇

\parSELECT

object_id

, ra

, dec

\parFROM

pdr2_wide.smallcat

LEFT JOIN pdr2_wide.specz USING (object_id)

LEFT JOIN pdr2_wide.photoz_mizuki USING (object_id)

WHERE

((photoz_best > 0.50 AND photoz_best <= 0.75) OR (specz_redshift > 0.50 AND
specz_redshift <= 0.75))

AND (i_kronflux_mag < 23 OR i_cmodel_mag < 23)

AND i_extendedness_value = 1

;

Figure 27: High-z Sample SQL query

##  Appendix C Additional Details About GaMPEN

To provide readers a visual understanding of how GaMPEN’s different
architectural components are organized, Figure 28 shows a schematic diagram
outlining the structure of both the STN and CNN in GaMPEN. For complete
details on individual layers in GaMPEN, please refer to Ghosh et al. (2022).

Figure 28: A schematic diagram of the Galaxy Morphology Posterior Estimation
Network. GaMPEN’s architecture consists of a downstream CNN module preceded by
an upstream STN module. The CNN module empowers GaMPEN to estimate posterior
distributions of galaxy morphology parameters. The upstream STN module trains
without any extra supervision and learns to apply appropriate cropping
transformations to the input image before passing it on to the CNN (for more
details about these modules, see §3). The numbers below each layer refer to
the number of filters/neurons in each layer. The yellow boxes inside the
convolutional layers show the kernel and the number beside it refers to the
corresponding kernel size. Only one kernel is shown per set of convolutional
layers; all other layers in the set have kernels of the same size. Conv2D and
ReLU refer to Convolutional Layers and Rectified Linear Units, respectively.

##  Appendix D Additional Details About Trained GaMPEN Models

As noted in §4, the training procedure of GaMPEN involves the tuning of
various hyper-parameters (e.g., learning rate, batch size, etc.). These hyper-
parameters are chosen based on the combination of the values that result in
the best performance on the validation data set. The final chosen hyper-
parameters for both the simulation trained GaMPEN models, as well as those
fine-tuned on real data, are shown in Table 4. Please refer to §4 for more
details on how we train these models.

Table 4: Tuned Values of Various Hyper-parameters  Model Name | Learning Rate | Momentum | L2 Regularization ($\lambda$) | Batch Size | Dropout Rate
---|---|---|---|---|---
Low-z Sim. Trained | $5\times 10^{-7}$ | 0.99 | $10^{-4}$ | 16 | $7\times 10^{-4}$
Mid-z Sim. Trained | $5\times 10^{-7}$ | 0.99 | $10^{-4}$ | 16 | $7\times 10^{-4}$
High-z Sim. Trained | $5\times 10^{-7}$ | 0.99 | $10^{-4}$ | 16 | $4\times 10^{-4}$
Low-z Final | $5\times 10^{-8}$ | 0.99 | $10^{-4}$ | 16 | $4\times 10^{-4}$
Mid-z Final | $5\times 10^{-8}$ | 0.99 | $10^{-4}$ | 16 | $2\times 10^{-4}$
High-z Final | $5\times 10^{-6}$ | 0.99 | $10^{-4}$ | 16 | $2\times 10^{-4}$

##  Appendix E Identifying Issues with our Light Profile Fitting Pipeline

We described in §5 a semi-automated pipeline that we used to determine the
structural parameters for $\sim 60,000$ galaxies using light-profile fitting.
After performing this analysis, we visually inspected the fits of a randomly
selected sub-sample of $\sim 300$ galaxies from each redshift bin, to assess
the quality of the fit.

The process of visually inspecting $\sim 900$ galaxies helped us to identify
three failure modes of our fitting pipeline:-

  1. 1.

for some galaxies, GALFIT assigned an extremely small axis ratio to one of the
components leading to an unphysical lopsided bulge/disk

  2. 2.

for some galaxies, the centroids of the two fitted components were too far
away from each other

  3. 3.

some galaxies had residuals that were too large

Note that for some galaxies, multiple failure modes were applicable. Examples
of each of these failure modes are shown in Figure 29. We use a combination of
different cuts on the fitted dataset to get rid of these failure modes – see
§5 for an extended discussion.

Figure 29: The different failure modes of the semi-automated light profile
fitting code described in §5. From left to right, we show the input image, the
mask generated by Source Extractor, the model generated by GALFIT, and the
residuals.

##  Appendix F Additional Two-Dimensional Residual & Uncertainty Plots

Figures 30 and 31 show the distribution of residuals (difference between
GaMPEN and GALFIT predictions) for the mid- and high-z bins.

Figures 32 and 33 show the uncertainties predicted by GaMPEN for each
parameter plotted against the predicted values for the mid- and high-z bins.

Figure 30: Residuals of the output parameters (difference between GaMPEN and
GALFIT predictions) plotted against the values predicted by GaMPEN for all
galaxies in the mid-z test set. See § 6.3 for details. Figure 31: Residuals of
the output parameters (difference between GaMPEN and GALFIT predictions)
plotted against the values predicted by GaMPEN for all galaxies in the high-z
test set. See § 6.3 for details. Figure 32: Uncertainties predicted by GaMPEN
for each parameter plotted against the predicted values for the mid-z test
set. The $\sigma$ for each parameter is defined as the width of the $68.27\%$
confidence interval. The line-shaded region in the top-left panel shows the
region where we recommend transforming quantitative $L_{B}/L_{T}$ predictions
to qualitative labels. See §6.4 for details. Figure 33: Uncertainties
predicted by GaMPEN for each parameter plotted against the predicted values
for the high-z test set. The $\sigma$ for each parameter is defined as the
width of the $68.27\%$ confidence interval. The line-shaded region in the top-
left panel shows the region where we recommend transforming quantitative
$L_{B}/L_{T}$ predictions to qualitative labels. See §6.4 for details.

## References

  * Aihara et al. (2018) Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, Publications of the Astronomical Society of Japan, 70, doi: 10.1093/pasj/psx066
  * Aihara et al. (2019) Aihara, H., Alsayyad, Y., Ando, M., et al. 2019, Publications of the Astronomical Society of Japan, 71, 114, doi: 10.1093/PASJ/PSZ103
  * Alam et al. (2015) Alam, S., Albareti, F. D., Prieto, C. A., et al. 2015, The Astrophysical Journal Supplement Series, 219, 12, doi: 10.1088/0067-0049/219/1/12
  * Bender et al. (1992) Bender, R., Burstein, D., & Faber, S. M. 1992, The Astrophysical Journal, 399, 462, doi: 10.1086/171940
  * Bernardi et al. (2013) Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, Monthly Notices of the Royal Astronomical Society, 436, 697, doi: 10.1093/mnras/stt1607
  * Bertin et al. (2020) Bertin, E., Schefer, M., Apostolakos, N., et al. 2020, in Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter (San Francisco: Astronomical Society of the Pacific), 461
  * Bertinl (1996) Bertinl, E. 1996, Astronomy and Astrophysics Supplement Series, 117, 393, doi: 10.1051/aas:1996164
  * Binney & Merrifield (1998) Binney, J., & Merrifield, M. 1998, Galactic astronomy (Princeton University Press), 796.  https://press.princeton.edu/titles/6358.html
  * Birrer & Amara (2018) Birrer, S., & Amara, A. 2018, Physics of the Dark Universe, 22, 189, doi: 10.1016/j.dark.2018.11.002
  * Blum & Mitchell (1998) Blum, A., & Mitchell, T. 1998, in Proceedings of the eleventh annual conference on Computational learning theory - COLT’ 98 (New York, New York, USA: ACM Press), 92–100, doi: 10.1145/279943.279962
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, Publications of the Astronomical Society of Japan, 70, doi: 10.1093/pasj/psx080
  * Bradshaw et al. (2013) Bradshaw, E. J., Almaini, O., Hartley, W. G., et al. 2013, Monthly Notices of the Royal Astronomical Society, 433, 194, doi: 10.1093/mnras/stt715
  * Chen et al. (2020) Chen, A., Chow, A., Davidson, A., et al. 2020, in Proceedings of the Fourth International Workshop on Data Management for End-to-End Machine Learning (New York, NY, USA: ACM), 1–4, doi: 10.1145/3399579.3399867
  * Cheng et al. (2021) Cheng, T.-Y., Conselice, C. J., Aragón-Salamanca, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 4425, doi: 10.1093/mnras/stab2142
  * Cool et al. (2013) Cool, R. J., Moustakas, J., Blanton, M. R., et al. 2013, The Astrophysical Journal, 767, 118, doi: 10.1088/0004-637X/767/2/118
  * Dekel et al. (2009a) Dekel, A., Sari, R., & Ceverino, D. 2009a, The Astrophysical Journal, 703, 785, doi: 10.1088/0004-637X/703/1/785
  * Dekel et al. (2009b) Dekel, A., Birnboim, Y., Engel, G., et al. 2009b, Nature 2009 457:7228, 457, 451, doi: 10.1038/nature07648
  * Dieleman et al. (2015) Dieleman, S., Willett, K. W., & Dambre, J. 2015, Monthly Notices of the Royal Astronomical Society, 450, 1441, doi: 10.1093/mnras/stv632
  * Dimauro et al. (2022) Dimauro, P., Daddi, E., Shankar, F., et al. 2022, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/MNRAS/STAC884
  * Dodge (2008) Dodge, Y. 2008, in The Concise Encyclopedia of Statistics (New York, NY: Springer New York), 502–505, doi: 10.1007/978-0-387-32833-1_379
  * Drinkwater et al. (2010) Drinkwater, M. J., Jurek, R. J., Blake, C., et al. 2010, Monthly Notices of the Royal Astronomical Society, 401, 1429, doi: 10.1111/j.1365-2966.2009.15754.x
  * Euclid Collaboration et al. (2022) Euclid Collaboration, Bretonnière, H., Kuchner, U., et al. 2022.  http://arxiv.org/abs/2209.12907
  * Ferrari et al. (2015) Ferrari, F., Carvalho, R. R. d., & Trevisan, M. 2015, The Astrophysical Journal, 814, 55, doi: 10.1088/0004-637X/814/1/55
  * Fomin et al. (2020) Fomin, V., Anmol, J., Desroziers, S., Kriss, J., & Tejani, A. 2020, High-level library to help with training neural networks in PyTorch, https://github.com/pytorch/ignite, GitHub
  * Garilli et al. (2014) Garilli, B., Guzzo, L., Scodeggio, M., et al. 2014, Astronomy & Astrophysics, 562, A23, doi: 10.1051/0004-6361/201322790
  * Genzel et al. (2008) Genzel, R., Burkert, A., Bouché, N., et al. 2008, The Astrophysical Journal, 687, 59, doi: 10.1086/591840
  * Ghosh et al. (2020) Ghosh, A., Urry, C. M., Wang, Z., et al. 2020, The Astrophysical Journal, 895, 112, doi: 10.3847/1538-4357/ab8a47
  * Ghosh et al. (2022) Ghosh, A., Urry, C. M., Rau, A., et al. 2022, The Astrophysical Journal, 935, 138, doi: 10.3847/1538-4357/ac7f9e
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature 2020 585:7825, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hausen & Robertson (2020) Hausen, R., & Robertson, B. E. 2020, The Astrophysical Journal Supplement Series, 248, 20, doi: 10.3847/1538-4365/ab8868
  * Haussler et al. (2007) Haussler, B., McIntosh, D. H., Barden, M., et al. 2007, The Astrophysical Journal Supplement Series, 172, 615, doi: 10.1086/518836
  * Häußler et al. (2022) Häußler, B., Vika, M., Bamford, S. P., et al. 2022, Astronomy & Astrophysics, 664, A92, doi: 10.1051/0004-6361/202142935
  * Hoyle (2016) Hoyle, B. 2016, Astronomy and Computing, 16, 34, doi: 10.1016/j.ascom.2016.03.006
  * Hsieh & Yee (2014) Hsieh, B. C., & Yee, H. K. C. 2014, The Astrophysical Journal, 792, 102, doi: 10.1088/0004-637X/792/2/102
  * Huertas-Company et al. (2015) Huertas-Company, M., Gravet, R., Cabrera-Vives, G., et al. 2015, The Astrophysical Journal Supplement Series, 221, 8, doi: 10.1088/0067-0049/221/1/8
  * Huertas-Company et al. (2016) Huertas-Company, M., Bernardi, M., Pérez-González, P. G., et al. 2016, Monthly Notices of the Royal Astronomical Society, 462, 4495, doi: 10.1093/mnras/stw1866
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ivezić et al. (2019) Ivezić, Z., Kahn, S. M., Tyson, J. A., et al. 2019, The Astrophysical Journal, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Kauffmann et al. (2004) Kauffmann, G., White, S. D. M., Heckman, T. M., et al. 2004, Monthly Notices of the Royal Astronomical Society, 353, 713, doi: 10.1111/j.1365-2966.2004.08117.x
  * Kawanomoto et al. (2018) Kawanomoto, S., Uraguchi, F., Komiyama, Y., et al. 2018, Publications of the Astronomical Society of Japan, 70, 66, doi: 10.1093/PASJ/PSY056
  * Kawinwanichakij et al. (2021) Kawinwanichakij, L., Silverman, J. D., Ding, X., et al. 2021, The Astrophysical Journal, 921, 38, doi: 10.3847/1538-4357/ac1f21
  * Kim & Brunner (2017) Kim, E. J., & Brunner, R. J. 2017, Monthly Notices of the Royal Astronomical Society, 464, 4463, doi: 10.1093/mnras/stw2672
  * Kloek & van Dijk (1978) Kloek, T., & van Dijk, H. K. 1978, Econometrica, 46, 1, doi: 10.2307/1913641
  * Kormendy (1979) Kormendy, J. 1979, The Astrophysical Journal, 227, 714, doi: 10.1086/156782
  * Kormendy & Kennicutt (2004) Kormendy, J., & Kennicutt, R. C. 2004, Annual Review of Astronomy and Astrophysics, 42, 603, doi: 10.1146/annurev.astro.42.053102.134024
  * Le Fèvre et al. (2013) Le Fèvre, O., Cassata, P., Cucciati, O., et al. 2013, Astronomy & Astrophysics, 559, A14, doi: 10.1051/0004-6361/201322179
  * Lilly et al. (2009) Lilly, S. J., Le Brun, V., Maier, C., et al. 2009, The Astrophysical Journal Supplement Series, 184, 218, doi: 10.1088/0067-0049/184/2/218
  * Liske et al. (2015) Liske, J., Baldry, I. K., Driver, S. P., et al. 2015, Monthly Notices of the Royal Astronomical Society, 452, 2087, doi: 10.1093/mnras/stv1436
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, 56–61, doi: 10.25080/Majora-92bf1922-00a
  * Momcheva et al. (2016) Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, The Astrophysical Journal Supplement Series, 225, 27, doi: 10.3847/0067-0049/225/2/27
  * Natarajan et al. (2021) Natarajan, P., Tang, K. S., Khochfar, S., et al. 2021, arXiv preprint arXiv:2103.13932, doi: 10.48550/arxiv.2103.13932
  * Newman et al. (2013) Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, The Astrophysical Journal Supplement Series, 208, 5, doi: 10.1088/0067-0049/208/1/5
  * Nishizawa et al. (2020) Nishizawa, A. J., Hsieh, B.-C., Tanaka, M., & Takata, T. 2020, Publ. Astron. Soc. Japan, 1, doi: 10.48550/arxiv.2003.01511
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, in Advances in Neural Information Processing Systems 32, ed. H. Wallach, H. Larochelle, A. Beygelzimer, F. d Alché-Buc, E. Fox, & R. Garnett (Curran Associates, Inc.), 8024–8035
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, The Astronomical Journal, 124, 266, doi: 10.1086/340952
  * Pentericci et al. (2018) Pentericci, L., McLure, R. J., Garilli, B., et al. 2018, Astronomy & Astrophysics, 616, A174, doi: 10.1051/0004-6361/201833047
  * Planck Collaboration et al. (2018) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, Astronomy and Astrophysics, 641, doi: 10.1051/0004-6361/201833910
  * Powell et al. (2017) Powell, M. C., Urry, C. M., Cardamone, C. N., et al. 2017, The Astrophysical Journal, 835, 22, doi: 10.3847/1538-4357/835/1/22
  * Pozzetti et al. (2010) Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, Astronomy & Astrophysics, 523, A13, doi: 10.1051/0004-6361/200913020
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, The Astronomical Journal, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Quionero-Candela et al. (2009) Quionero-Candela, J., Sugiyama, M., Schwaighofer, A., & Lawrence, N. D. 2009, Dataset Shift in Machine Learning (The MIT Press)
  * Racca et al. (2016) Racca, G. D., Laureijs, R., Stagnaro, L., et al. 2016, in Proceedings of SPIE Astronomical Telescopes + Instrumentation, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, 99040O, doi: 10.1117/12.2230762
  * Recht et al. (2019) Recht, B., Roelofs, R., Schmidt, L., & Shankar, V. 2019.  http://arxiv.org/abs/1902.10811
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy & Astrophysics, 558, A33, doi: 10.1051/0004-6361/201322068
  * Robotham et al. (2017) Robotham, A. S. G., Taranu, D. S., Tobar, R., Moffett, A., & Driver, S. P. 2017, Monthly Notices of the Royal Astronomical Society, 466, 1513, doi: 10.1093/mnras/stw3039
  * Rowe et al. (2015) Rowe, B. T., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121, doi: 10.1016/j.ascom.2015.02.002
  * Schawinski et al. (2007) Schawinski, K., Kaviraj, S., Khochfar, S., et al. 2007, The Astrophysical Journal Supplement Series, 173, 512, doi: 10.1086/516631
  * Schawinski et al. (2014) Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, Monthly Notices of the Royal Astronomical Society, 440, 889, doi: 10.1093/mnras/stu327
  * Sellwood (2014) Sellwood, J. A. 2014, Reviews of Modern Physics, 86, 1, doi: 10.1103/RevModPhys.86.1
  * Shallue & Vanderburg (2018) Shallue, C. J., & Vanderburg, A. 2018, The Astronomical Journal, 155, 94, doi: 10.3847/1538-3881/aa9e09
  * Sharma et al. (2020) Sharma, K., Kembhavi, A., Kembhavi, A., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 2280, doi: 10.1093/MNRAS/STZ3100
  * Shimakawa et al. (2021) Shimakawa, R., Tanaka, T. S., Toshikage, S., & Tanaka, M. 2021, Publications of the Astronomical Society of Japan, 73, 1575, doi: 10.1093/pasj/psab097
  * Silverman et al. (2015) Silverman, J. D., Kashino, D., Sanders, D., et al. 2015, The Astrophysical Journal Supplement Series, 220, 12, doi: 10.1088/0067-0049/220/1/12
  * Simard et al. (2011) Simard, L., Trevor Mendel, J., Patton, D. R., Ellison, S. L., & McConnachie, A. W. 2011, Astrophysical Journal, Supplement Series, 196, 11, doi: 10.1088/0067-0049/196/1/11
  * Simard et al. (2002) Simard, L., Willmer, C. N. A., Vogt, N. P., et al. 2002, The Astrophysical Journal Supplement Series, 142, 1, doi: 10.1086/341399/FULLTEXT/
  * Simmons & Urry (2008) Simmons, B. D., & Urry, C. M. 2008, The Astrophysical Journal, 683, 644, doi: 10.1086/589827
  * Spergel et al. (2013) Spergel, D., Gehrels, N., Breckinridge, J., et al. 2013.  http://arxiv.org/abs/1305.5422
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., & Salakhutdinov, R. 2014, Journal of Machine Learning Research, 15, 1929, doi: 10.5555/2627435
  * Tanaka (2015) Tanaka, M. 2015, The Astrophysical Journal, 801, 20, doi: 10.1088/0004-637X/801/1/20
  * Tanaka et al. (2022) Tanaka, T. S., Shimakawa, R., Shimasaku, K., et al. 2022, Publications of the Astronomical Society of Japan, 74, 1, doi: 10.1093/pasj/psab105
  * Taori et al. (2020) Taori, R., Dave, A., Shankar, V., et al. 2020, in Advances in Neural Information Processing Systems, ed. H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, & H. Lin, Vol. 33 (Curran Associates, Inc.), 18583–18599.  https://proceedings.neurips.cc/paper/2020/file/d8330f857a17c53d217014ee776bfd50-Paper.pdf
  * Tarsitano et al. (2022) Tarsitano, F., Bruderer, C., Schawinski, K., & Hartley, W. G. 2022, Monthly Notices of the Royal Astronomical Society, 511, 3330, doi: 10.1093/mnras/stac233
  * Tarsitano et al. (2018) Tarsitano, F., Hartley, W. G., Amara, A., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 2018, doi: 10.1093/mnras/sty1970
  * Tremaine et al. (2002) Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, The Astrophysical Journal, 574, 740, doi: 10.1086/341002
  * Tuccillo et al. (2018) Tuccillo, D., Huertas-Company, M., Decencière, E., et al. 2018, Monthly Notices of the Royal Astronomical Society, 475, 894, doi: 10.1093/mnras/stx3186
  * van der Wel (2008) van der Wel, A. 2008, The Astrophysical Journal, 675, L13, doi: 10.1086/529432
  * Van Der Wel et al. (2012) Van Der Wel, A., Bell, E. F., Häussler, B., et al. 2012, Astrophysical Journal, Supplement Series, 203, 24, doi: 10.1088/0067-0049/203/2/24
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, The Astrophysical Journal, 788, 28, doi: 10.1088/0004-637X/788/1/28
  * Vega-Ferrero et al. (2021) Vega-Ferrero, J., Domínguez Sánchez, H., Bernardi, M., et al. 2021, Monthly Notices of the Royal Astronomical Society, 506, 1927, doi: 10.1093/mnras/stab594
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walmsley et al. (2020) Walmsley, M., Smith, L., Lintott, C., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 1554, doi: 10.1093/mnras/stz2816
  * Walmsley et al. (2021) Walmsley, M., Lintott, C., Géron, T., et al. 2021, Monthly Notices of the Royal Astronomical Society, 509, 3966, doi: 10.1093/mnras/stab2093
  * Weinmann et al. (2006) Weinmann, S. M., van den Bosch, F. C., Yang, X., & Mo, H. J. 2006, Monthly Notices of the Royal Astronomical Society, 366, 2, doi: 10.1111/j.1365-2966.2005.09865.x
  * Wuyts et al. (2011) Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, The Astrophysical Journal, 742, 96, doi: 10.1088/0004-637X/742/2/96
  * York et al. (2000) York, D. G., Adelman, J., Anderson, Jr., J. E., et al. 2000, The Astronomical Journal, 120, 1579, doi: 10.1086/301513
  * Zanisi et al. (2021) Zanisi, L., Huertas-Company, M., Lanusse, F., et al. 2021, Monthly Notices of the Royal Astronomical Society, 501, 4359, doi: 10.1093/mnras/staa3864